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CO ' Abstract 

00 ' 

In this paper we describe Herwig++ version 2.3, a general-purpose Monte Carlo event 
generator for the simulation of hard lepton-lepton, lepton-hadron and hadron-hadron col- 



lisions. A number of important hard scattering processes are available, together with an 
qq ■ interface via the Les Houches Accord to specialized matrix element generators for addi- 

, tional processes. The simulation of Beyond the Standard Model (BSM) physics includes a 

range of models and allows new models to be added by encoding the Feynman rules of the 
model. The parton-shower approach is used to simulate initial- and final-state QCD radia- 
tion, including colour coherence effects, with special emphasis on the correct description of 
radiation from heavy particles. The underlying event is simulated using an eikonal multiple 
parton-parton scattering model. The formation of hadrons from the quarks and gluons 
produced in the parton shower is described using the cluster hadronization model. Hadron 
decays are simulated using matrix elements, where possible including spin correlations and 
off-shell effects. 
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1 Introduction 



Herwig++ is a general-purpose event generator for the simulation of high-energy lepton-lepton, 
lepton-hadron and hadron-hadron collisions with special emphasis on the accurate simulation of 
QCD radiation. It builds upon the heritage of the HERWIG program [1-6], while providing a 
much more flexible structure for further development. It already includes several features more 
advanced than the last FORTRAN version. Herwig++ provides a full simulation of high energy 
collisions with the following special features: 

• Initial- and final-state QCD jet evolution taking account of soft gluon interference via 
angular ordering; 

• A detailed treatment of the suppression of QCD radiation from massive particles, the dead- 
cone effect [7]; 

• The simulation of BSM physics including correlations between the production and decay 
of the BSM particles together with the ability to add new models by simply encoding the 
Feynman rules; 

• An eikonal model for multiple partonic scatterings to describe the underlying event [8]; 

• A cluster model of the hadronization of jets based on non-perturbative gluon splitting; 

• A sophisticated model of hadron and tau decays using matrix elements to give the momenta 
of the decay products for many modes and including a detailed treatment of off-shell effects 
and spin correlations. 

Some of these features were already present in the first version of Herwig++ [9]. However, there 
have been many improvements to both the physics and structure of the simulation following 
this first release, most notably the extension to hadron-hadron collisions. Given the significant 
differences between the current version of the program, 2.3, and that described in [9] we will 
describe all of the features of the program in this paper. 

A number of other generators are also being (re-)written for the LHC era. The PYTHIA 
event generator is currently being rewritten as PYTHIA8 [10]. The rewrite of ARIADNE [11] is in 
progress as well. Like Herwig++, this is built on the platform of ThePEG [12], which we describe 
below. SHERPA [13] is a completely new event generator project. 

It is useful to start by recalling the main features of a generic hard, high-momentum transfer, 
process in the way it is simulated by Herwig++. The processes involved can be divided into a 
number of stages corresponding to increasing time and distance scales: 

1. Elementary hard subprocess. In the hard process the incoming particles interact to pro- 
duce the primary outgoing fundamental particles. This interaction can involve either the 
incoming fundamental particles in lepton collisions or partons extracted from a hadron in 
hadron-initiated processes. In general this is computed at leading order in perturbation 
theory, although some processes calculated at next-to-leading order are now included [14] 
and work is ongoing to include additional processes [15-17]. The energy scale of the hard 
process, together with the colour flow between the particles, sets the initial conditions for 
the production of QCD radiation in the initial- and final-state parton showers. 
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2. Initial- and final-state parton showers. The coloured particles in the event are perturba- 
tively evolved from the hard scale of the collision to the infrared cutoff. This occurs for 
both the particles produced in the collision, the final-state shower, and the initial partons 
involved in the collision for processes with incoming hadrons, the initial-state shower. The 
coherence of the emission of soft gluons in the parton showers from the particles in the hard 
collision is controlled by the colour flow of the hard collision. Inside the parton shower, it is 
simulated by the angular ordering of successive emissions. The choice of evolution variable 
together with the use of quasi-collinear splitting functions allows us to evolve down to zero 
transverse momentum for the emission, giving an improved simulation of the dead-cone 
effect for radiation from massive particles [7]. 

3. Decay of heavy objects. Massive fundamental particles such as the top quark, electroweak 
gauge bosons, Higgs bosons, and particles in many models of physics beyond the Standard 
Model, decay on time-scales that are either shorter than, or comparable to that of the QCD 
parton shower. Depending on the nature of the particles and whether or not strongly inter- 
acting particles are produced in the decay, these particles may also initiate parton showers 
both before and after their decay. One of the major features of the Herwig++ shower 
algorithm is the treatment of radiation from such heavy objects in both their production 
and decay. Spin correlations between the production and decay of such particles are also 
correctly treated. 

4. Multiple scattering. For large centre-of-mass energies the parton densities are probed in a 
kinematic regime where the probability of having multiple partonic scatterings in the same 
hadronic collision becomes significant. For these energies, multiple scattering is the domi- 
nant component of the underlying event that accompanies the main hard scattering. These 
additional scatterings take place in the perturbative regime, above the infrared cut-off, and 
therefore give rise to additional parton showers. We use an eikonal multiple scattering 
model [8], which is based on the same physics as the FORTRAN JIMMY package [18], to- 
gether with some minor improvements. In addition to that we included non-perturbative 
partonic scatters below the infrared cut-off [19], which enables us to simulate minimum 
bias events as well as the underlying event in hard scattering processes. 

5. Hadronization. After the parton showers have evolved all partons involved in hard scatter- 
ings, additional scatters and partonic decays down to low scales, the final state typically 
consists of coloured partons that are close in momentum space to partons with which they 
share a colour index, their colour 'partner' (in the large N c limit this assignment is unique). 
Herwig++ uses the cluster hadronization model [2] to project these colour-anticolour pairs 
onto singlet states called clusters, which decay to hadrons and hadron resonances. The 
original model of Ref. [2], which described this decay as pure phase space has been progres- 
sively refined as described in Sect. 7. Clusters that are too massive or too light for decay 
directly to hadrons to provide a good description are treated differently, again described in 
Sect. 7. 

6. Hadron decays. The hadron decays in Herwig++ are simulated using a matrix element 
description of the distributions of the decay products, together with spin correlations be- 
tween the different decays, wherever possible. The treatment of spin correlations is fully 
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integrated with that used in perturbative production and decay processes so that corre- 
lations between the production and decay of particles like the tau lepton, which can be 
produced perturbatively but decays hadronically, can be treated consistently. 

The program and additional documentation are available from 

http : // pro j ects . hepf orge . org/herwig 

This manual concentrates on the physics included in the Herwig++ simulation, which has been 
the subject of a number of publications [9,14-17,20-26]. Additional documentation of the code, 
together with examples of how to use the program and further information is available from our 
website and wiki. We provide a bug-tracker, which should be used to report any problems with 
the program or to request user support. 

Herwig++ is distributed under the GNU General Public License (GPL) version 2. This 
ensures that the source code will be available to users, grants them the freedom to use and 
modify the program and sets out the conditions under which it can be redistributed. However, it 
was developed as part of an academic research project and is the result of many years of work by 
the authors, which raises various issues that are not covered by the legal framework of the GPL. 
It is therefore distributed together with a set of guidelines 1 , agreed by the MCnet collaboration, 
which set out various expectations that we have of responsible users. In particular, concerning 
citation of relevant physics publications, they state that the main software reference as designated 
by the program authors (i.e. this manual for Herwig++ versions 2.1 onwards) should always be 
cited, as well as the original literature on which the program is based to the extent that it is of 
relevance for a study, applying the same threshold criteria as for other literature. To help users 
in this, Herwig++ produces a ETgX file that lists the primary physics citation(s) for each module 
that has been active during a given run. The authors are always happy to help users determine 
which citations are relevant to a particular study. 

The remainder of this manual is set out as follows. The next section contains a brief technical 
description which should be sufficient to understand the details of the program included in the 
discussion of the physics simulation. More detailed technical documentation can be obtained 
from the website above, including Doxygen descriptions of all classes. 

The rest of the manual then discusses the physics of each stage of the simulation process in 
detail, describing the physics models used in the simulation, together with the main parameters 
of the models and the structure of the code. Finally, we give a summary and our plans for 
future improvements. Appendices give some more technical information, a series of examples of 
the program in use, and a brief description of the process by which the default parameters were 
tuned to data. 



These guidelines are contained in the GUIDELINES file distributed with the release and are also available from 
http : //www.montecarlonet . org/ index .php?p=Publications/Guidelines 
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2 Technical Details 



While this manual is primarily a description of the physics models used in Herwig++, by its 
nature we cannot wholly avoid discussing the technical details of the program. We need to 
discuss some aspects of the program's structure and the mechanism for changing physics model 
parameters, so that users can adjust parameters, change the hard process they are simulating, 
or make any of the other modifications that are necessary to make the program useful to an 
individual user. In this section we will give a basic overview of the structure of the program, 
which is designed to supplement the Doxygen documentation of the source code available at 

http : //projects .hepforge . org/herwig/doxygen 

Herwig++ is based on ThePEG [12] — the Toolkit for High Energy Physics Event Generation, 
a framework for implementing Monte Carlo event generators. ThePEG provides all parts of the 
event generator infrastructure that do not depend on the physics models used as a collection of 
modular building blocks. The specific physics models of Herwig++ are implemented on top of 
these. 

Each part of Herwig++ is implemented as a C++ class that contains the implementation of 
the Herwig++ physics models, inheriting from an abstract base class in ThePEG. This allows the 
implementations of different physics models to live side-by-side and be easily exchanged. 

The central concept in ThePEG is the Repository, which holds building blocks in the form of 
C++ objects that can be combined to construct an EventGenerator object, which in turn will be 
responsible for all steps of event generation. Within the Repository, one can create objects, set 
up references between them, and change all parameter values. The Repository object needs to be 
populated with references to all required objects for the physics models used at run time. The 
objects can then be persistently stored, or combined to produce an EventGenerator. The default 
Repository layout for Herwig++ is shown in Table 1. The composition of the Repository is con- 
trolled through a simple configuration language, described in Appendix A. This set of commands 
allows the user to configure the generator at run time. Through this mechanism, selection of 
different physics models or different model parameters is possible without recompilation of the 
code. 

The EventGenerator object is responsible for the run 2 as a whole. It holds the infrastructure 
objects that are needed for the run, like the generation of random numbers, the particle properties 
stored as ParticleData objects, and handles any exceptions. 

The actual generation of each event is the responsibility of the EventHandler. It manages the 
generation of the hard scattering process 3 and the subsequent evolution of the event through five 
StepHandler objects, each of which is responsible for generating one main part of the event: 

1. The SubProcessHandler is responsible for generating the hard sub-process as described in 
Section 3. This handler is skipped if the hard process is read in from a Les Houches Accord 
event file. 

2. The CascadeHandler generates the parton shower from the hard process. 
2 The generation of a series of events. 

3 The generation of the hard process by the EventHandler and its inheriting classes is discussed in more detail 
in Section 3.5. 
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Table 1: Overview of the default Repository layout for Herwig++. Each box represents a reference 
to an independent C++ object held in the repository, which can be swapped out for a different 
implementation. 



3. The MultiplelnteractionHandler produces additional hard scatters when using a multiple 
parton-parton scattering model to simulate the underlying event in hadron-hadron colli- 
sions. In practice, given the close relationship between the parton shower and the additional 
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hard scatters in Herwig++, the multiple scattering model is implemented as part of the 
Herwig++ implementation of the CascadeHandler, the ShowerHandler. 



4. The Had ronization Handler is responsible for the formation of hadrons out of the quarks and 
gluons left after the parton shower. 

5. The DecayHandler is responsible for decaying both the unstable hadrons produced by the 
HadronizationHandler, and any unstable fundamental particles that may have been produced 
in either the hard process or parton shower. 

The StepHandler base classes in ThePEG do not implement any physics models themselves. 
This must be done by inheriting classes, which provide an implementation of a specific model. 
The Herwig++ ShowerHandler for example, inherits from CascadeHandler and implements the 
Herwig++ parton shower model by overriding the virtual cascadeQ member function. 

In addition to the five main handlers, the EventHandler allows for pre- and post-handlers to 
be called before and after each step. This allows for additional processing of the event where 
required: in Herwig++ BSM physics or top quark production, the HwDecayHandler is used as 
a pre-handler for the ShowerHandler to ensure that all the unstable fundamental particles have 
decayed before the parton shower occurs. 

The implementation of a physics model as a StepHandler generally does not put all the code 
needed for the simulation in one class, but makes use of an, often large, number of helper classes. 

This brief description only discusses the classes responsible for generating the core parts of 
the event. Other classes and concepts are discussed in more detail in the Doxygen documentation. 

The mechanisms for exploring and changing the values of switches and parameters are also 
described in Appendix A. It is worth mentioning that 'default' values of switches and parameters 
can appear in one of two places: the repository entries in the default . in files; or the class 
constructors and at present there is no built-in mechanism to ensure that they are consistent. 
When both exist, the former takes precedence. The values described as 'default' in this manual 
are those that appear in the default . in files. A further confusion appears, because the value 
described as default in the Doxygen documentation is not guaranteed to be the same as either of 
the others. A mechanism to ensure that all three default values are the same will be introduced in 
a future version, but until then, users are reminded that the default . in files remain the primary 
source of parameter values. 
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3 Matrix Elements 



In Herwig++ the library of matrix elements for QCD and electroweak processes is relatively small, 
certainly with respect to the large range of processes available in its FORTRAN predecessor [5,6]. 
Indeed, the library of Standard Model processes is largely intended to provide a core of important 
processes with which to test the program. Whereas, at the time of the development of the 
original FORTRAN program, matrix elements needed to be calculated and implemented by hand, 
nowadays there are a number of programs that automate these calculations, for a wide range of 
processes with high multiplicity final states. It has therefore been our intention that, in general, 
users should study most processes of interest via our interface to these programs. 

Nevertheless, there are still some cases for which it is useful to have Herwig++ handle all 
stages of the event generation process. This is particularly true for processes in which spin 
correlations between the production and decay stages are significant e.g. those involving top 
quarks or tau leptons. Such correlation effects are hard to treat correctly if different programs 
handle different steps of the simulation process. 

In order to facilitate the process of adding new matrix elements, where needed, and to enable 
us to generate the spin correlation effects [27-30] , we have based all matrix element calculations 
on the helicity libraries of ThePEG. As well as providing a native library of Standard Model 
processes and an interface to parton-level generators, Herwig++ also includes matrix elements 
for hard 2^2 collisions and 1 — > 2 and 1 — >• 3 decays, arising in various models of new physics 
(see Sect. 5). 

Starting with version 2.3 a number of next-to-leading-order (NLO) matrix elements in the 
POsitive Weight Hardest Emission Generator (POWHEG) scheme of Refs. [31, 32] are also in- 
cluded. 

3.1 Leading order matrix elements 

For e + e~ colliders only four hard processes are included: 

• Quark- ant iquark production, via interfering photon and Z° bosons, is implemented in the 
MEee2gZ2qq class. No approximation is made regarding the masses of the particles. This 
process is essential for us to validate the program using QCD analyses of LEP data. 

• Dilepton pair production, via interfering photon and Z° bosons, is implemented in the 
MEee2gZ2ll class. No approximation is made regarding the masses of the particles 4 . This 
process is used to check the implementation of spin correlations in r decays. 

• The Bjorken process, Z°h° production, which is implemented in the MEee2ZH class. This 
process is included as it is very similar to the production of Z°h° and W ± h° in hadron- 
hadron collisions and uses the same base class for most of the calculation. 

• The vector-boson fusion (VBF) processes, e + e~ — > e + e~h° and e + e~ — > u e i? e h , are imple- 
mented in the MEee2HiggsVBF class. 

For deep inelastic scattering (DIS) only two processes are included. Neutral and charged 
current processes are implemented in the MENeutralCurrentDIS and MEChargedCurrentDIS classes, 

4 i-channel photon and Z° boson exchange are not included. 
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respectively. In neutral current processes both the incoming and outgoing partons are considered 
to be massless, whereas in the charged current process the masses of the outgoing partons are 
included. For neutral current scattering both photon and Z° boson exchange are included. 

A much wider range of matrix elements is included in the standalone code for the simulation 
of events in hadron colliders: 

• Difermion production via s-channel electroweak gauge bosons. The matrix elements for the 
production of fermion-antifermion pairs through W ± bosons, or interfering photons and Z° 
bosons, are implemented in the MEqq2W2ff and MEqq2gZ2ff classes respectively. Only 
s-channel electroweak gauge boson diagrams are included for the hadronic modes. 

• The production of a Z° or W boson in association with a hard jet is simulated using 
the MEPP2ZJet or MEPP2WJet class respectively. The decay products of the bosons are 
included in the 2 — > 3 matrix element and the option of including the photon for Z° 
production is supported. 

• The 2^2 QCD scattering processes are implemented in the MEQCD2to2 class. Currently 
all the particles are treated as massless in these processes. 

• The matrix element for the production of a heavy quark- ant iquark pair (top or bottom 
quark pairs), is coded in the MEPP2QQ class. No approximations are made regarding the 
masses of the outgoing qq pair. 

• The MEPP2GammaGamma class implements the matrix element for the production of 
prompt photon pairs. In addition to the tree-level qq — > 77 process the loop-mediated 
99 ~^ 77 process is included. 

• Direct photon production in association with a jet is simulated using the MEPP2GammaJet 
class. As with the QCD 2 — > 2 process all of the particles are treated as massless in these 
processes. 

• The production of an s-channel Higgs boson via both gg — > h° and qq — > h° is simulated 
using the MEPP2Higgs class. 

• The production of a Higgs boson in association with the Z° or W ± bosons is simulated 
using the MEPP2ZH or MEPP2WH class respectively. 

• The production of the Higgs boson in association with a hard jet is simulated using the 
MEPP2HiggsJet class. 

In addition we have a matrix element class, MEQCD2to2Fast, that uses hard-coded formulae 
for the QCD 2^2 scattering matrix elements rather than the helicity libraries of ThePEG. This 
class is significantly faster than the default MEQCD2to2 class, although it does not implement 
spin correlations. It is intended to be used in the generation of the multiple parton-parton 
scatterings for the underlying event where the spin correlations are not important but due to 
the number of additional scatterings that must be generated the speed of the calculation can 
significantly affect the run time of the event generator. There is also the MEMinBias class which 
is only used to simulate soft scattering processes as part of the underlying event model. 
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3.2 Next-to-leading-order matrix elements 

In recent years there have been a number of additional developments which aim to improve on the 
results of parton shower simulations by providing a description of the hardest emission together 
with a next-to-leading order (NLO) cross section [15-17,31-41] 5 . 

The first successful scheme for matching at NLO was the MCONLO approach [33-38] which 
has been implemented with the HERWIG event generator for many processes. The method has 
two draw backs; first, it involves the addition of correction terms that are not positive definite 
and therefore can result in events with a negative weight and second, the implementation of the 
method is heavily dependent on the details of the parton shower algorithm used by the event 
generator. 

In Ref. [31] a novel method, referred to as POWHEG (POsitive Weight Hardest Emission 
Generator), was introduced to achieve the same aims as MCONLO while creating only positive 
weight events and being independent of the event generator with which it is implemented. The 
POWHEG method has been applied to Z pair hadropro duct ion [39], heavy flavour hadroproduc- 
tion [41], e + e~ annihilation to hadrons [15] and top production in lepton collisions [17]. A general 
outline of the ingredients required for POWHEG with two popular NLO subtraction schemes is 
given in Ref. [32]. 

The POWHEG shower algorithm involves generating the hardest emission in separately 
using a Sudakov form factor containing the full matrix element for the emission of an extra 
parton and adding to this vetoed showers, which produce radiation at lower scales in the shower 
evolution variable, and a truncated shower, which generates radiation at higher scales in the 
shower evolution variable, than the scale of the highest pt emission. While the POWHEG scheme 
is independent of the parton shower algorithm, it does require the parton shower to be able to 
produce vetoed and truncated showers. The ability to perform vetoed showers is present in most 
modern Monte Carlo event generators, however some changes are required to enable them to 
generate truncated showers. Although the POWHEG approach is formally correct to the same 
accuracy as the M CO NLO technique the two methods differ in their treatment of sub-leading 
terms. 

A small number of processes [14] in the POWHEG scheme are now implemented in Herwig++ 
together with a full implementation of the truncated shower. These processes are implemented 
in the following way: 

• the matrix elements are calculated with NLO accuracy and a Born configuration supplied 
in the same way as for the leading-order matrix elements; 

• a modified PowhegEvolver is used to generate the shower which uses the relevant 
HardestEmissionGenerator to generate the hardest emission from the shower; 

• the event is then showered, including the truncated shower, as described in Sect. 6.9. 
Currently the following processes are implemented: 

• the Drell Yan production of neutral vector bosons 7/Z is simulated using the 
MEqq2gZ2ffPowheg class; 

5 There have been other theoretical ideas but only the MCONLO and POWHEG methods have led to practical 
programs whose results can be compared with experimental data. 
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• the Drell Yan production of charged vector bosons, i.e. W^, is implemented in the 
MEqq2W2ffPowheg class; 

• the production of the Higgs boson via the gluon-gluon fusion process is simulated using the 
MEPP2HiggsPowheg class; 

• the production of the Higgs boson in association with the boson is implemented in the 
MEPP2WHPowheg class; 

• the production of the Higgs boson in association with the Z° boson is simulated using the 
MEPP2ZHPowheg. 

The hardest emission for these processes is then generated using the DrellYanHardGenerator or 
GGtoHHardGenerator for vector boson and Higgs boson production respectively. More details of 
the simulation of QCD radiation can be found in Sect. 6.9 and Ref. [14]. 

3.3 Les Houches interface 

There are a number of matrix element generators available that can generate parton-level events 
using either the original Les Houches Accord [42] or the subsequent extension [43], which spec- 
ified a file format for the transfer of the information between the matrix element generator and 
a general-purpose event generator, such as Herwig++, rather than the original FORTRAN COM- 
MON block. 

In addition to the internal mechanism for the generation of hard processes, ThePEG provides 
a general LesHouchesEventHandler class, which generates the hard process using the Les Houches 
Accord. In principle a run-time interface could be used to directly transfer the information 
between the matrix element generator and Herwig++, however we expect that the majority of 
such interfaces will be via data files containing the event information using the format specified 
in Ref. [43]. 

We expect that this approach will be used for the majority of hard processes in Herwig++. 

3.4 Processes with incoming photons 

It is possible to have hard scattering processes with incoming photons in hadron-hadron collisions, 
for example in the higher-order QED corrections to the Drell- Yan production of W ± or Z° bosons. 
These can not be directly showered by the Herwig++ parton shower and therefore we provide 
a IncomingPhotonEvolver class which can be used as one of the PreCascadeHandlers in these 
processes to perform the backward evolution of the photon to a quark or antiquark which can 
then be evolved by the Herwig++ parton shower. 

This performs one backward branching evolving in transverse momentum from a starting 
scale PTstart given by the pt of the softest particle in the event, or a minimum scale minpT if 
the scale is below the minimum allowed value. This is performed using a Sudakov form factor, 

\ Jpl Pt 2 W* V l xf y (x,p' T ) 

where pt is the transverse momentum of the branching, a is the fine structure constant, x is 
the momentum fraction of the photon, e, is the electric charge of the particle produced in the 
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backward evolution and the sum over i runs over all the quarks and antiquarks. The splitting 
function is 



where z is the fraction of the momentum of incoming parton produced in the backward evolution 
given to the photon. The pr and momentum fraction of the branching are generated in the 
same way as those in the parton shower, as described in Sect. 6. The momenta of the particles, 
including the new branching are then reconstructed as described in Sect. 6.5. 

3.5 Code structure 

In ThePEG the generation of the hard process is the responsibility of the EventHandler. The base 
EventHandler class only provides the abstract interfaces for the generation of the hard process with 
the actual generation of the kinematics being the responsibility of inheriting classes. There are 
two such classes provided in ThePEG: the Standard EventHandler, which implements the internal 
mechanism of ThePEG for the generation of the hard process; and the LesHouchesEventHandler, 
which allows events to be read from data files. 

3.5.1 StandardEventHandler 

The StandardEventHandler uses a SubProcessHandler to generate the kinematics of the particles 
involved in the hard process. In turn the SubProcessHandler makes use of a number of MEBase 
objects to calculate the matrix element and generate the kinematics for specific processes. The 
specific matrix elements used in a given run of the EventGenerator can be specified using the 
MatrixElements interface of the SubProcessHandler. The MEBase object is responsible for: 

• defining the particles that interact in a given process, by specifying a number of DiagramBase 
objects; one DiagramBase is specified per flavour combination. 

• returning the differential partonic cross section 



when supplied with the partonic centre-of-mass energy of the collision and n random num- 
bers between and 1. Each MEBase class specifies how many random numbers it requires 
to calculate the partonic cross section and kinematics for the processes it implements. For 
example a 2 — > 2 process typically needs two 6 random numbers, one each for the polar and 
azimuthal angles. 

• creating a HardVertex object describing the interaction that occurred, including the spin- 
unaveraged matrix element to allow spin correlation effects to be generated. 

One MEBase object is generally used to describe one physical process with different partonic 
flavours. The selection of flavours within each subprocess is carried out internally by the 

6 In practice as the matrix elements do not depend on the azimuthal angle we often only use one random 
number for the polar angle and generate the second random number locally. 




(3.2) 



z 
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EventHandler. The resulting cross sections can be output with varying levels of detail, con- 
trolled by the StatLevel switch; by default they are only broken down by MEBase objects. The 
SubProcessHandler then uses a SamplerBase object to perform the unweighting of the cross section 
and generate events with unit weight. In practice for 2^2 cross section the generation of the 
kinematics and other technical steps is handled by the ME2to2Base class. In addition the actual 
calculation of the matrix element can be easily implemented using the Helicity classes of ThePEG. 
All of the matrix elements in Herwig++ inherit 7 from ME2to2Base and make extensive use of the 
Helicity library of ThePEG. 

In general the main switch for the generation of the hard process is the 
MatrixElements interface, which allows the MEBase objects to be specified and hence de- 
termines which hard scattering processes are generated. In addition, each class inheriting from 
MEBase in Herwig++ has a number of parameters that control the incoming, outgoing and inter- 
mediate particles in a specific process. These are controlled by Interfaces in the specific matrix 
element classes. A number of different partonic subprocesses can be handled at the same time 
by simply specifying several MEBase objects. 

3.5.2 LesHouchesEventHandler 

The LesHouchesEventHandler class inherits from the EventHandler class of ThePEG. The class has 
a list of LesHouchesReader objects that are normally connected to files with event data produced 
by an external matrix element generator program, although it could in principle include a direct 
run-time link to the matrix element generator or read events 'on the fly' from the output of a 
matrix element generator connected to a pipe. 

When an event is requested by LesHouchesEventHandler, one of the readers is chosen according 
to the cross section of the process for which events are supplied by that reader. An event is read 
in and subsequently handled in the same way as for an internally generated process. The use of 
the LesHouchesEventHandler class is described in Appendix B.10. 

3.5.3 Kinematic cuts 

For cuts on the hard process we use Cuts objects from ThePEG. All cuts applied to the generation 
of the hard process can be specified via its Interfaces. There are many types of cuts that can be 
applied. 

Cuts applied to the overall hard process, such as a minimum or maximum invariant mass M 
of the process, can be specified directly as a parameter of the Cuts class. The minimum value of 
the invariant mass for the hard process is set using the MhatMin parameter. Many more cuts 
can be specified by using the Interfaces of the Cuts class. Among those that are used in Herwig++ 
are cuts on the momentum fractions of the incoming partons and the hard process scale. 
The default set of cuts we apply in hadronic collisions is M > 20 GeV (MhatMin), i li2 > 1CT 5 
(XlMin, X2Min) and Q > 1 GeV (ScaleMin). 

In addition to these general cuts it is possible to specify cuts that are only applied to particular 
particles, particle pairs or resonant intermediate particles. In order to do so, one has to specify a 
number of OneCutBase, TwoCutBase or MultiCutBase objects in the Cuts object that is applied. 

7 The only exception is the MEQCD2to2Fast class, which is 'hand written' for speed. 
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Whenever we use OneCutBase cuts we use either the SimpleKTCut class for massless particles 
or the KTRapidityCut class for massive particles. These require that a Matcher object is set up for 
the particles to which the cut is applied. The Matcher classes used in Herwig++ all inherit from 
MatcherType. In addition to the Matcher classes provided by ThePEG, Herwig++ provides addi- 
tional matchers for top quarks (TopMatcher), photons (Photon Matcher), W ± (WBoson Matcher), 
Z° (ZBosonMatcher) and Higgs (HiggsBosonMatcher) bosons. This can be either a single parti- 
cle, for example the top quark, or a group of particles, like the leptons. Then, for example, the 
minimum transverse momentum of that particle /c_i_ im m can be specified as MinKT. In addition 
we use minimum and maximum values of pseudorapidity via MinEta and MaxEta for massless 
particles using the SimpleKTCut class or MinRapidity and MaxRapidity for massive particles 
using the KTRapidityCut class. By default we use k i ) TT1 i n > 20 GeV for partons and \rj\ < 3 for 
photons from the hard scattering process. 

An example of a MultiCutBase class is the V2LeptonsCut class. We use it to limit the invariant 
mass of lepton pairs. It is given similarly to the general cut as MinM. We use the rather loose 
cut 20 GeV < M < 1.4 TeV by default. Another useful parameter of this class is the specification 
of the lepton families (Families) or the charge combination (CComb) of the lepton pair the cut 
is applied to. 

As the cuts are applied to all the particles produced in the collision, for W^/Z production in 
association with either a jet or a Higgs boson the cuts are also applied to the decay products of 
the boson. This can lead to inefficiencies in the generation of the hard process and a suppression 
of the hadronic boson decays with the default cuts on the quarks. 
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4 Perturbative Decays and Spin Correlations 



In Herwig++ the decays of the fundamental particles and the unstable hadrons are handled 
in the same way in order that correlation effects for particles such as the tau lepton, which is 
produced perturbatively but decays non-perturbatively, are correctly treated. Eventually it is 
intended that the unstable fundamental particles will be decayed during the parton-shower stage 
of the event, however currently in order that the correlation effects are correctly generated all 
the perturbative particle decays are performed before the generation of the parton shower by 
using the HwDecayHandler as one of the PreCascadeHandlers in the EventHandler responsible 
for generating the event. The Decayer classes used in Herwig++ to perform the decays of the 
fundamental Standard Model particles make use of the Helicity classes of ThePEG to calculate 
the helicity amplitudes for the decay matrix elements. The code structure for the Decayer classes 
used in Herwig++ and the HwDecayHandler implement the algorithm of Refs. [27-30] to correctly 
include the spin correlations. 

In the next subsection we describe the spin correlation algorithm of [27-30] using the example 
of top production and decay. This is followed by a description of the modelling of the decay of the 
fundamental particles of the Standard Model, the production and decays of particles in models 
of physics Beyond the Standard Model is discussed in Section 5. We then describe the simulation 
of QED radiation in particle decays. Finally we briefly discuss the structure of the code for the 
decays of fundamental particles. 



4.1 Spin correlations 

When calculating the matrix element for a given hard process or decay one must take into account 
the effect of spin correlations, as they will affect the distributions of particles in the final state. 
In particular these correlations are important in the production and decay of the top quark, for 
the production of tau leptons in Higgs decays and in models of BSM physics where one can have 
two models that possess a very similar particle spectrum but with particles that have different 
spins. 

An algorithm for correctly incorporating these correlations into a Monte Carlo is demonstrated 
in Refs. [27-30]. Rather than discuss the algorithm in full detail here we will describe it by 
considering the example of the process e + e~ — > tt where the top quark subsequently decays, via 
a W + boson, to a b quark and a pair of light fermions. 

Initially, the outgoing momenta of the tt pair are generated according to the usual cross- 
section integral 

(2tt) 4 f d 3 p t d 3 P t +e _^ tt - e+e 



2s J (2^2E t {2nf2E- t Ml ^ ~* " 

where M.^- ^ s the ma trix element for the initial hard process and A 4 f are the helicities of 
the t and t respectively. One of the outgoing particles is then picked at random, say the top, and 
a spin density matrix calculated 



Pli>. = jfMZ£-**M$f**, (4.2) 



with iV defined such that Tr p = 1. 
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The top is decayed and the momenta of the decay products distributed according to 



(2tt) 4 f d 3 p b d 3 Pw+ , ut ^ bW+u ^ bW + , s 



2m t J (2ir) 3 2E b (2ir) 3 2E w + rx t x t A < A ^+ W+ 

where the inclusion of the spin density matrix ensures the correct correlation between the top 
decay products and the beam. 

A spin density matrix is calculated for the W + only, because the b is stable 

PT = Ip\ ,'M^M^ w+ , (4.4) 

and the W + decayed in the same manner as the top. Here the inclusion of the spin density matrix 
ensures the correct correlations between the W + decay products, the beam and the bottom quark. 

The decay products of the W + are stable fermions so the decay chain terminates here and a 
decay matrix for the W + 

DT. V = ^MSSMt?**, (4.5) 



is calculated. Moving back up the chain a decay matrix for the top quark is calculated using the 
decay matrix of the W + , 

D\ x , = ^ T Mt$ w +M?r? w+ DY + w • (4.6) 

A t A t N A t A w -+ x t X w+ X w+ X w+ V I 

Since the top came from the hard scattering process we must now deal with the t in a similar man- 
ner but instead of using 5\ t \' when calculating the initial spin density matrix, the decay matrix of 
the top is used and the i decay is generated accordingly. The density matrices pass information 
from one decay chain to the associated chain thereby preserving the correct correlations. 

The production and decay of the top, using the spin correlation algorithm, is demonstrated in 
Figs. 1-3. The hard scattering process and subsequent decays were generated using the general 
matrix elements described in Sect. 5 rather than the default ones. The results from the full matrix 
element calculation are also included to show that the algorithm has been correctly implemented. 
The separate plots illustrate the different stages of the algorithm at work. Figure 1 gives the 
angle between the beam and the outgoing lepton. The results from the simulation agree well 
with the full matrix element calculation, which demonstrates the consistency of the algorithm 
for the decay of the t. 

Figure 2 gives the angle between the top quark and the produced lepton. This shows the same 
agreement as the previous figure and demonstrates the correct implementation of the spin den- 
sity matrix for the t decay. Finally, Fig. 3 gives the results for the angle between the final-state 
lepton/anti- lepton pair showing the correct implementation of the decay matrix that encodes the 
information about the t decay. Distributions for various processes within the Minimal Super- 
symmetric Standard Model and for tau production in Higgs decay are shown in Refs. [24,25]. 

The same algorithm is used regardless of how the particles are produced, in order to consis- 
tently implement the spin correlations in all stages of the event generation process. 
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Figure 1: Angle between the beam and the outgoing lepton in e + e~ — > ti — > bbl + vil~tii in the lab frame 
for a centre-of-mass energy of 500 GeV with (a) unpolarized incoming beams, (b) negatively polarized 
electrons and positively polarized positrons and (c) positively polarized electrons and negatively polar- 
ized electrons. The data points show the results of the simulation as production and decay including 
spin correlations, while the histograms use the full matrix element for e + e _ — > ti — > bbl + vil~Vi. 
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Figure 2: Angle between the lepton and the top quark in e + e~ — > ti — > bbl + vil~vi in the lab frame for 
a centre-of-mass energy of 500 GeV with (a) unpolarized incoming beams, (b) negatively polarized elec- 
trons and positively polarized positrons and (c) positively polarized electrons and negatively polarized 
electrons. The data points show the results of the simulation as production and decay including spin 
correlations, while the histograms use the full matrix element for e + e~ — > ti —> bbl + vil~i>i. 
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Figure 3: Angle between the outgoing lepton and anti-lepton in e + e~ — ► ti — ► bbl^ v\l~ v\ in the lab 
frame for a centre-of-mass energy of 500 GeV with (a) unpolarized incoming beams, (b) negatively 
polarized electrons and positively polarized positrons and (c) positively polarized electrons and nega- 
tively polarized electrons. The data points show the results of the simulation as production and decay 
including spin correlations, while the histograms use the full matrix element for e + e~ — ► ti — ► bbl + vil~Ui. 

4.2 Standard Model decays 

There are a small number of decays of fundamental Standard Model particles currently imple- 
mented. These are implemented as Decayer classes for top quark, W ± and Z°, and Higgs boson 
decays. The following classes are available: 

• the SMTopDecayer implements the three-body decay of the top quark to the bottom quark 
and a Standard Model fermion-antifermion pair, via an intermediate W + boson; 

• the SMWZ Decayer class implements the decay of the W ± and Z° bosons to a Standard 
Model fermion-antifermion pair; 

• the SMHiggsFermionsDecayer class implements the decay of the Higgs boson to a Standard 
Model fermion-antifermion pair, i.e. h° — > //; 

• the SMHiggsWWDecayer implements the decay of the Higgs boson to W ± or Z° bosons, 
i.e. h° —> W + W~, Z°Z°, including the decay of the gauge bosons; 

• the SMHiggsGGHiggsPPDecayer implements the decay of the Higgs boson to a pair of either 
gluons or photons. 

In many cases off-shell effects for the electroweak gauge bosons are included by generating 
the gauge bosons as intermediate particles, for example in top quark and Higgs boson decays. 
In general, external top quarks and W and Z° bosons are produced off mass-shell using the 
approach described in Ref. [44]. The Higgs boson mass is generated in the same way as in the 
FORTRAN HERWIG program using the more sophisticated approach described in Ref. [45]. 
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4.3 QED radiation 



The simulation of QED radiation using the approach of Ref. [23] has been included for both 
particle decays and unstable s-channel resonances produced in the hard process. This approach 
is based on the YFS formalism [46] , which takes into account large double- and single- soft photon 
logarithms to all orders. In addition, the leading collinear logarithms are included to O (a) by 
using the dipole splitting functions. By default the production of QED radiation is switched 
off for both decays and hard processes. It may be included by using the QEDRadiationHandler 
in the EventHandler as one of the PostSubProcessHandlers for the hard process or using the 
PhotonGenerator interface of the relevant Decayer inheriting from the Decay Integrator class for 
the decays. 

4.4 Code structure 

The code structure for particle decays in Herwig++ is described in more detail in Sect. 9.6 
for the hadronic decays. All of the Decayer classes for fundamental particles inherit from the 
Decayl integrator class in order to use the multi-channel phase-space integration it provides. 

The SMHiggsMassGenerator implements the generation of the mass of off-shell Higgs bosons 
using the running width implemented in the SMHiggsWidthGenerator class. These classes inherit 
from the GenericMassGenerator and GenericWidthGenerator classes of Herwig++ in order to have 
access to the full infrastructure for the simulation of off-shell particles described in Sect. 9. 

The structure of the code for the simulation of QED radiation in particle decays is designed 
to be general, so that other approaches can be implemented. The generation of the radiation 
is handled by a class inheriting from the abstract DecayRadiationGenerator class. Currently only 
the YFS approach, as described in Ref. [23], is implemented in the SOPHTY class, which uses the 
helper FFDipole and IFDipole classes for radiation from final-final and initial-final dipoles, respec- 
tively. In addition the QEDRadiationHandler is included to allow the DecayRadiationGenerator to 
be used to generate radiation in the decay of particles generated as s-channel resonances in the 
hard process. 
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5 Physics Beyond The Standard Model 



No one knows what kind of physics will be encountered in the LHC era and it is likely that a 
variety of new physics models will need to be considered in determining its exact nature. This 
eventuality has been accounted for in the design of the Herwig++ program, by the inclusion of 
a general framework for the implementation of new physics models. Using this framework, new 
models can be realized quickly and efficiently This method is described in full in Refs. [24, 44] 
and will be reviewed here. 

In describing the features needed to simulate Beyond the Standard Model (BSM) processes, we 
need only concern ourselves with the hard collisions, either producing known particles through 
modified couplings or the exchange of new particles, or producing new particles in the final 
state, and with decays of the new particles. All other steps of event generation are handled 
in the same way as for Standard Model processes 8 . Both of these steps involve calculating an 
amplitude, which in turn relies on knowledge of the Feynman rules within the model being used. 
In Herwig++ the Feynman rules are implemented as a series of Vertex classes, which inherit from 
the generic classes of ThePEG. These Vertex classes are based on the HELAS formalism [47], with 
each class able to evaluate the vertex as a complex number or, given different information, an 
off-shell wavefunction that can be used as input for another calculation. Each Feynman diagram 
contributing to a given process is evaluated in terms of these vertex building blocks and the sum 
of the resulting contributions is squared to give the matrix element. 

In this section we start by briefly describing the generation of the hard processes and decays 
in models of new physics, this is followed by a description of models currently implemented in 
Herwig++, including the Standard Model, and the structure of the code. 

5.1 Hard process 

Section 3 gave details on the default matrix elements available for generating Standard Model 
processes in Herwig++. These classes are based on specific particle interactions whereas the 
classes used for BSM models are based on the external spin structure of a 2 — > 2 scattering 
process. To generate a specific process the user specifies the desired states that are to participate 
in the hard interaction, using the configuration files, and the code then generates the relevant 
diagrams and a MatrixElement object for each process 9 . 

The generic matrix elements use a colour flow decomposition to calculate the value of |.M| 2 . 
This method cuts down on the amount of colour algebra necessary in the evaluation of QCD 
processes by rewriting the colour structures of certain diagrams in terms of others in the same 
process. As an example, consider the process q a — > g c g d , which has diagrams with amplitudes 
given by 



8 Othcr features do emerge in certain models, for example the hadronization of new long-lived coloured particles, 
which is not yet fully implemented in HerwigH — h, but for the majority of new physics models under active study 
this is the case. 

9 It is only necessary to specify a single outgoing particle as the code will produce all processes with this particle 
in the final state. 




(5.1a) 
(5.1b) 
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where Ai{t, u ,s} is the colour-stripped amplitude for each diagram type. Using the colour matrix 
identities, Eq. (5.1c) can be rewritten as [t c ,t d ]b a Ai s and is then a combination of the other two 
colour structures. By defining a colour flow /, as a combination of colour-stripped amplitudes 
possessing the same colour structure, in this case f\ = M. t — -Ms and f'2 = M. u + -Ms, we 
can cut down the number of colour factors that need to be evaluated. The full matrix element 
squared, summed over final-state spins and colours and averaged over initial spins and colours, 
is obtained by adding up products of colour flows and the appropriate colour factor. For any 
process ab —> cd this can be written as 

\M\ 2 = Z^^Y, C ^ff ( 5 - 2 ) 
b a b b C a C b ^ 

where Cjj is a matrix containing the squared colour factors, denotes the ith colour flow for 
the set of helicities A, Z is an identical particle factor, S a ,b is the number of polarization states 
for each incoming particle and C a ^ is the number of colour states for each incoming particle. 

To carry out the parton showering and hadronization stages of the simulation we must assign 
a colour to each particle participating in each hard collision. This information is needed in 
determining the initial conditions for the parton shower (Sect. 6.3), and how clusters are formed 
in the hadronization model (Sect. 7). To this end, each fundamental coloured particle is associated 
to a ColourLine object. For the particles involved in the hard interactions, the colour assignments 
are made by selecting a colour flow from a list contained in the corresponding MatrixElement class 
as follows. Once a momentum configuration for the primary hard scattering has been generated, 
each colour flow is assigned a weight according to how much it contributes to the total value of 
the matrix element (neglecting the interference between them, which is typically suppressed by 
l/iV 2 and also by dynamical effects). One of these colour flows is then probabilistically chosen 
on the basis of this weight distribution. 



5.2 Decays 

To be able to decay the BSM states, the possible decay modes must first be known. If a su- 
persymmetric model is required one can use a spectrum generator to produce not only the 
required spectrum, in accordance with the SUSY Les Houches Accord [48], but also a decay ta- 
ble. Herwig++ is designed to be able to read this information and set up the appropriate decay 
modes for later use. Other models do not have such programs and therefore the list of possible 
two- and three-body decays is generated automatically. 

When generating the possible decays automatically we also need to be able to calculate the 
partial width of a given mode so that the branching fraction and total width can be calculated. 
For a general two-body decay, the matrix element only depends on the mass-square values of each 
particle so the phase-space factor can be integrated separately and the partial width is given by 



r (a -> b, c) = — — , 5.3 



where \Ai\ 2 is the matrix element squared summed over final-state colours and spins and averaged 
over initial-state colours and spins and p cm is the centre-of-mass momentum 

Pcm = ^-[H- (m 6 + m c ) 2 ) (ml - (m b - m c ) 2 )] 1/2 . (5.4) 
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A three-body decay has a partial width given by 



1 
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(5.5) 



a J (m b +m c ) 2 



where 
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(5.6a) 
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(5.6b) 



with E* = {m 2 ab — m 2 + m 2 c )/2mb c and E* d = (m 2 — m\ c — rrij)/2mb c . In general, the phase-space 
integration can no longer be performed analytically since the matrix element is a complicated 
function of the invariant mass combinations rribc and vn c ^ therefore it must be performed numer- 
ically. Given the low number of dimensions of the phase-space integrals in Eq. (5.5), they are 
performed using standard techniques rather than by the Monte Carlo method. The total width 
of the parent is simply the sum of the partial widths. 

To compute the momenta of the decay products we need to be able to calculate the matrix 
element for a selected decay mode. When each mode is created it is assigned a Decayer object 
that is capable of calculating the value of \Ai\ 2 for that process. It is done in a similar manner 
to the hard matrix element calculations, i.e. using the helicity libraries of ThePEG. 

In decays involving coloured particles that have more than one possible colour flow, the colour 
is treated in exactly the same way as described in Sect. 5.1 for hard processes. 



The production and decay processes described above have their external particles on mass shell 
throughout. This assumes that the narrow width approximation, defined by the following as- 



1. the resonance has a small width T compared with its pole mass M, T <C M; 

2. we are far from threshold, — M ^> T, where ^/s denotes the centre-of-mass energy; 

3. the propagator is separable; 

4. the mass of the parent is much greater than the mass of the decay products; 

5. there are no significant non-resonant contributions; 

is a valid approximation. In general, given that we do not have a specific mass spectrum, this is 
not a good enough approximation. In particular if processes occur at or close to threshold, there 
can be large corrections that we need to take into account. 

To improve our simulation we provide an option to include the weight factor 



5.3 Off-Shell Effects 



sumptions: 




(5.7) 
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throughout the production and decay stages, where T(m) is the running width of the particle 
to be considered off shell, M is the pole mass and m min , max are defined such that the maximum 
deviation from the pole mass is a constant times the on-shell width. A derivation of this factor 
can be found in the appendix of Ref . [44] . 

5.4 Model descriptions 

This section will give a description of the models that are included in Herwig++. In general in 
Herwig++ the implementation of a physics model consists of a main class, which inherits from 
the Standard Model class and implements the calculation of any parameters required by the model 
or, for a SUSY model, reads them from an input SUSY Les Houches file. In addition, there 
are various classes that inherit from the general Vertex classes of ThePEG, which implement the 
Feynman rules of the model. There may also be some classes implementing other features of the 
model, for example the running couplings in the specific model. 

5.4.1 Standard Model 

The implementation of the Standard Model in Herwig++ inherits from the 
StandardModelBase class of ThePEG. ThePEG includes classes to implement the running strong 
and electromagnetic couplings, together with the CKM matrix. 

In Herwig++ we include our own implementations of the running electromagnetic coupling, 
in the AlphaEM class, and the running strong coupling in the 02AlphaS class. By default we use 
the implementations of the running couplings from ThePEG and the Herwig++ implementations 
are only provided to allow us to make exact comparisons with the FORTRAN HERWIG program. 

In order to perform helicity amplitude calculations we need access to the full CKM matrix. 
However the CKM Base class of ThePEG only provides the squares of the matrix elements. The 
StandardCKM class therefore provides access to the matrix elements as well and it is used in all 
our helicity amplitude calculations. 

We have also included a structure for the implementation of running mass calculations. The 
RunningMassBase class provides a base class and the two-loop QCD running mass is implemented 
in the RunningMass class. 

The Standard Model input parameters in Herwig++ do not form a minimal set in that it is 
possible to independently set the value of the weak mixing angle in such a way that the tree-level 
relationship between the W ± and Z° boson masses is not satisfied. The electroweak parameters 
we use are: 

coupling at zero momentum transfer, 
=0.232]; 

and Z°, Mz = 91.1876 GeV, bosons, which are 

• the mixing angles, 12 [theta_12=0.2262], 13 [theta_13=0.0037] and 
6> 23 [theta_23=0.0413], and phase, S [delta=1.05], of the CKM matrix. 

In addition, many of the Standard Model couplings to the Z° boson can be changed to simulate 
non-Standard Model effects if desired. 



• the value of the electromagnetic 
[EW/AlphaEM=137.04]; 

• the value of sm 2 6 w , [EW/Sin2ThetaW= 

• the masses of the W ± , M w = 80.403 GeV, 
taken from their ParticleData objects: 
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Spin 


Particles 





d L ,u L ,s Ll c L ,bi,ti 

dR, Ur, Sr, Cr, b 2 ,t 2 

H°, A , H+ 


1/2 


9i X?) X25 X*]) X45 Xi 1 X2 



Table 2: The additional particle content of the MSSM contained in Herwig++. The particle's 
PDG codes are the standard ones given by the Particle Data Group [49]. 



5.4.2 Minimal Supersymmetric Standard Model 

The Minimal Supersymmetric Standard Model (MSSM) is the most studied supersymmetric 
model and as such it should be included in any generator attempting to simulate BSM physics. 
As its name suggests it contains the smallest number of additional fields required for the theory 
to be consistent. The additional particle content over that of the Standard Model is listed in 
Table 2. 

The additional particles must have masses and couplings to be of any use in an event simu- 
lation. For supersymmetric models various programs are available that, given some set of input 
parameters, produce a spectrum containing all of the other parameters necessary to be able to 
calculate physical quantities within the model. As stated in the previous section the output from 
such a generator must comply with the SUSY Les Houches Accord (SLHA) [48] for it to be used 
with Herwig++. 

While reading the information from an SLHA file is straightforward, there is a minor compli- 
cation when dealing with particle masses that have a mixing matrix associated with them. For 
example, consider the neutralinos, which are an admixture of the bino b, third wino W3 and 2 
higgsinos hi and h 2 . The physical eigenstates Xi are given by 

Xi = Nitfl (5.8) 

where Nij is the neutralino mixing matrix in the ip° = (—ib, —iw, hi, h 2 ) T basis. The diagonalized 
mass term for the gauginos is then N*ftAr N\ which in general can produce complex mass values. 
To keep the mass values real the phase is instead absorbed into the definition of the corresponding 
field thereby yielding a strictly real mass and mixing matrix. There is however a price to be paid 
for this — while the masses are kept real they can become negative. For an event generator a 
negative mass for a physical particle does not make sense so we instead choose a complex-valued 
mixing matrix along with real and non-negative masses. If a negative mass is encountered while 
reading a Les Houches file, the physical mass is taken as the absolute value and the appropriate 
row of the mixing matrix is multiplied by a factor of i. This approach is used in order to facilitate 
the implementation of extended supersymmetric models in the future. 

5.4.3 Randall-Sundrum Model 

The first models proposed with extra dimensions were of the Randall-Sundrum (RS) [50] type 
where a tensor particle, namely the graviton, is included and is allowed to propagate in the extra 
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dimensions. All other matter, however, is restricted to our usual 4D brane and as a result all of 
the SM couplings are left unchanged. The only extra couplings required are those of the graviton 
to ordinary matter, which depend on a single parameter A„.. 

5.4.4 Minimal Universal Extra Dimensions Model 

We also include a model based on the idea of universal extra dimensions where all fields are 
allowed to propagate in the bulk. Following similar lines to supersymmetry, the model included 
in Herwig++ is of a minimal type and has a single compact extra dimension of radius R [51]. 

Compactifying the extra dimension and allowing all fields to propagate in it leads to a rich new 
structure within the theory. Analogous to the particle-in-a-box scenario, one obtains an infinite 
number of excitations of the fields all characterized by a quantity called the KK-number. This is 
most easily demonstrated by showing how a scalar field $ would decompose after compactification 



where are the 4D coordinates, y the position in the 5th dimension and n is the KK-number of 
the mode with n = identified as the SM mode. In general, in some compactification schemes, 
it is possible to have KK-number-violating interactions but in the Minimal Universal Extra 
Dimensions (MUED) framework in Herwig++ we include only those interactions that conserve 
KK-parity P = (—1)™ and also limit ourselves to n = 1. 

Table 3 shows the MUED particle content contained in Herwig++ along with their particle 
ID codes as these have not been standardized by the Particle Data Group [49]. Unlike the MSSM 
there are no external programs available that calculate the mass spectrum so this must be done 
internally by the UEDBase class, which inherits from the StandardModel class and implements the 
UED model. At tree level the mass of any level-n particle is simply given by (m^ + (n/R) 2 ) 1 ^ 2 , 
where m is the mass of the SM particle, and 1/R is generally much larger than the SM mass so 
the spectrum is highly degenerate and no decays can occur. This situation changes once radiative 
corrections are taken into account and a spectrum that can be phenomenologically similar to the 
MSSM arises. The full set of radiative corrections, as derived in Ref. [52], is incorporated in the 
UEDBase class to give a realistic spectrum. 

5.5 Code structure 

The ModelGenerator class is responsible for setting up the new MatrixElement objects, which 
inherit from the GeneralHardME class, and DecayMode objects for a new physics model. Helper 
classes aid in the creation of these objects, they are: 

HardProcessConstructor the HardProcessConstructor is responsible for creating the dia- 
grams for the requested processes and constructing the appropriate GeneralHardME object(s); 

ResonantProcessConstructor the ResonantProcessConstructor is of a similar design to the 
HardProcessConstructor but it only constructs the resonant diagrams for a process; 




(5.9) 
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Particle 


ID code 


Spin 


Particle 


ID code 





h\ 
Al 

Ht 


5100025 
5100036 
5100037 


1 


9*i 

7i* 
z o* 

W+* 


5100021 
5100022 
5100023 
5100024 


1/2 


dl 
ul 

s'i 

bl 
t'l 


5100001 
5100002 
5100003 
5100004 
5100005 
5100006 


1/2 


d° ± 
ul 
si 

bl 
tl 


6100001 
6100002 
6100003 
6100004 
6100005 
6100006 




eV 
Ki 
A*i ' 

— • 

Ki 


5100011 
5100012 
5100013 
5100014 
5100015 
5100016 




— 

— o 

r l 


6100011 
6100013 
6100015 



Table 3: The MUED particle spectrum contained in Herwig++ along with their ID codes. * 
denotes a doublet under SU(2) and ° a singlet. As with the standard PDG codes an antiparticle 
is given by the negative of the number in the table. 



DecayConstructor the DecayConstructor stores a collection of objects that inherit from the 
NBodyDecayConstructor class. Each of these is responsible for constructing the decay modes 
for the n-body decays. Currently the TwoBodyDecayConstructor class, for two-body decays, the 
ThreeBodyDecayConstructor class for three-body decays and the WeakCurrentDecayConstructor 
class, for weak decays using the weak currents from Sect. 9.3.1 for decays where two particles are 
almost mass degenerate, are implemented. 

In addition, the ModelGenerator class is responsible for setting up objects of 
BSMWidthGenerator and MassGenerator type so that off-shell effects can be simulated. To achieve 
this either, ParticleData objects are added to the Offshell interface so that the selected particles 
are treated as off shell, or the WhichOffshell interface is set to All so that all BSM particles 
are treated as off shell. 

The matrix element classes all inherit from the Genera IHard ME class and implement the matrix 
element for a particular spin configuration. The classes inheriting from the GeneralHardME class 
and the spin structures they implement are given in Table 4. 

The on-shell decayer classes inherit from either the Genera ITwoBodyDecayer or 
GeneralThreeBodyDecayer class and each is responsible for calculating the value of the matrix 
element for that particular set of spins. A GeneralTwoBodyCurrentDecayer class also exists for de- 
cay modes created with the WeakCurrentDecayConstructor class. The Decayer classes implemented 
in Herwig++ and the types of decay they implement are given in Table 5. 

The use of BSM physics models is described in Appendix B.8 where examples of using all the 
models included with the release are given. 

The specification of the particles involved in the hard process is achieved through the 
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Class Name 


Hard Process 


MEfTCff 
MEff2ss 
MEff2vs 
MEff2vv 
MEfv2fs 
MEfv2vf 
MEvv2ff 
MEvv2ss 
MEvv2vv 


Fermion fermion to fermion fermion. 
Fermion fermion to scalar scalar. 
Fermion fermion to vector scalar. 
Fermion fermion to vector vector. 
Fermion vector to fermion scalar. 
Fermion vector to vector fermion. 
Vector vector to fermion fermion. 
Vector vector to scalar scalar. 
Vector vector to vector vector. 



Table 4: The general hard process matrix elements, based on spin structures, implemented in 
Herwig++. 



Class Name 


Decay 


FFSDecayer 


Fermion to fermion scalar decay. 


FFVDecayer 


Fermion to fermion vector decay. 


FFVCurrentDecayer 


Fermion to fermion vector decay with the vector off-shell 
and decaying via a weak current from Sect. 9.3.1. 


FtoFFFDecayer 


Fermion to three fermion decay. 


FtoFVVDecayer 


Fermion to fermion and two vector decay. 


SFFDecayer 


Scalar to fermion fermion decay. 


SSSDecayer 


Scalar to two scalar decay. 


SSVDecayer 


Scalar to scalar vector decay. 


SVVDecayer 


Scalar to two vector decay. 


SVVLoopDecayer 


Scalar to two vector decay via a loop. 


StoSFFDecayer 


Scalar to scalar and two fermion decay. 


StoFFVDecayer 


Scalar to two fermion and vector decay. 


VFFDecayer 


Vector to two fermion decay. 


VSSDecayer 


Vector to two scalar decay. 


VVVDecayer 


Vector to two vector decay. 


VtoFFVDecayer 


Vector to two fermion and vector decay. 


TFFDecayer 


Tensor to two fermion decay. 


TSSDecayer 


Tensor to two scalar decay 


TVVDecayer 


Tensor to two vector decay. 



Table 5: The general decays based on spin structures implemented in Herwig+- K 
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Class 


Interaction 


SMFFGVertex 
SMFFPVertex 
SMFFWVertex 
SMFFZVertex 
SMFFHVertex 
SMGGGVertex 
SMGGGGVertex 
SMWWWVertex 
SMWWWWVertex 
SMWWHVertex 
SMHGGVertex 
SMHPPVertex 


Interaction of the gluon with the SM fermions 
Interaction of the photon with the SM fermions 
Interaction of the boson with the SM fermions 
Interaction of the Z° boson with the SM fermions 
Interaction of the Higgs boson with the SM fermions 
Triple gluon vertex 
Four gluon vertex 

Triple electroweak gauge boson vertex 
Four electroweak gauge boson vertex 

Interaction of the Higgs boson with the electroweak gauge bosons 

Higgs boson coupling to two gluons via quark loops 

Higgs boson coupling to two photons via fermion and boson loops 



Table 6: Herwig++ Vertex classes for the Standard Model. 

Incoming and Outgoing interfaces of the HardProcessConstructor. Both interfaces are lists 
of ParticleData objects. The switch IncludeEW can be set to No to include only the strong 
coupling diagrams. 

In order to pass spin correlations through the decay stage, Decayl integrator objects must 
be created. This is achieved by populating a list held in the ModelGenerator class, which can 
be accessed through the DecayParticles interface. The particles in this list will have spin 
correlation information passed along when their decays are generated. If a decay table is read in 
for a SUSY model then the CreateDecayModes interface should be set to No so that only the 
decay modes listed in the externally generated decay table are created 10 . For all other models 
the possible decay modes are also created from the particles in the DecayParticles list. 

In addition to the code that handles the calculation of the matrix elements for the decays 
and scattering cross sections each model requires a number of classes to implement the model. 

The Standard Model is implemented in the Standard Model class, which inherits from the 
StandardModelBase class of ThePEG and implements access to the helicity Vertex classes and 
some additional couplings, such as the running mass, used by Herwig++. The Vertex classes that 
implement the Standard Model interactions are given in Table 6. 

The structure of the implementation of the MSSM in Herwig++ is designed to allow extended 
SUSY models to be added in the future. Therefore the SusyBase class, which inherits from the 
Standard Model class, is designed to read in the SLHA files specifying the SUSY spectrum. The 
details of the MSSM are implemented in the MSSM class, which inherits from the SusyBase class. 
The Vertex classes for the MSSM are given in Table 7. A spectrum file in SLHA format must be 
supplied, as described in Appendix B.8.1, or the MSSM model cannot be used. 

The RSModel class inherits from the StandardModel class and implements the calculations 
needed for the Randall-Sundrum model. We have only implemented the vertices that are phe- 
nomenologically relevant and therefore some four-point vertices that are not important for res- 
onance graviton production are not included. The Vertex classes implemented for the Randall- 
Sundrum model are given in Table 8. 

10 If a decay table is being used with a SUSY model then the DecayParticles list must still be populated so that 
the decays will have spin correlation information included. 
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Class 


Interaction 


SS NFS Vertex 


Neutralino with a SM fermion and a sfermion 


SSCFSVertex 


Chargino with a SM fermion and a sfermion 


SSGFSVertex 


Gluino with a quark and squark 


SSNNZVertex 


A pair of neutralinos with a boson 


SSCCZVertex 


A pair of charginos with a Z° boson 


SSCNWVertex 


Chargino with a neutralino and a boson 


SSGSGSGVertex 


SM gluon with a pair of gluinos 


SSGSSVertex 


SM gluon with a pair of squarks 


SSWSSVertex 


SM gauge boson with a pair of sfermions 


SSFFHVertex 


A pair of SM fermions with a Higgs boson 


SSWHH Vertex 


SM electroweak gauge bosons with a pair of Higgs bosons 


SSWWHVertex 


A pair of gauge bosons with a Higgs boson 


SSGOGOHVertex 


A pair of gauginos with a Higgs boson 

! O O OO 


SSHSFSFVertex 


A Higgs boson with a pair of sfermions 


SSHHH Vertex 


Triple Higgs boson self coupling 


SSHGGVertex 


A Higgs boson with a pair of gluons via quark and squark loops 


SSGGSQSQVertex 


A pair of gluons with a pair of squarks 



Table 7: Herwig++ Vertex classes for the MSSM. 



Class 


Interaction 


RSModelFFGRVertex 


Coupling 


of the 


graviton to 


SM 


fermions 


RSModelSSGRVertex 


Coupling 


of the 


graviton to 


the 


Higgs boson 


RSModelFFVGRVertex 


Coupling 


of the 


graviton to 


two 


SM 




fermions 


and a £ 


jauge boson 




RSModelVVGRVertex 


Coupling 


of the 


graviton to 


two 


gauge bosons 


RSModelVVVGRVertex 


Coupling 


of the 


graviton to 


three gauge bosons 



Table 8: Herwig++ Vertex classes for the Randall-Sundrum model. 

Two parameters can be controlled in the Randall-Sundrum model; the cutoff A,,, and the 
mass of the graviton. The default mass of the graviton is 500 GeV and this can be changed via 
the NominalMass interface of its ParticleData object. The cutoff is set via the Lambda_pi 
interface of the RSModel object and has a default value of lOTeV. 

The UED model is implemented in the UEDBase class, which inherits from the Standard Model 
class and implements the calculation of the parameters of the model. The Vertex classes for the 
UED model are given in Table 9. 

There are three parameters that can be set to control the UED model: the inverse of the radius 
of compactification i? -1 ; the cutoff scale A; and the mass of the Higgs boson at the boundary of 
the compactified dimension m^. These are controlled through the interfaces: 

InverseRadius the value of R~ l , the default value is 500 GeV; 

LambdaR the dimensionless number AR, the default value is 20; 
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Class 


Interaction 


UEDFlFlPOVertex 
UEDFlFlWOVertex 
UEDFlFlZOVertex 
UEDFlFlGOVertex 
UEDFlFOWlVertex 
UEDFlFOGlVertex 
UEDFlFOHlVertex 
UEDPOHlHlVertex 
UEDWOWlWlVertex 
UEDGlGlGOVertex 
UEDGOGOGlGlVertex 
UEDWOA1 HI Vertex 

UEDZOHlHlVertex 
UEDZOAlhlVertex 


SM photon with a pair of KK-1 fermions 

SM boson with a pair of KK-1 fermions 

SM Z° boson with a pair of KK-1 fermions 

SM gluon with a pair of KK-1 fermions 

KK-1 fermion with an EW KK-1 boson and a SM fermion 

KK-1 fermion with a KK-1 gluon and a SM fermion 

KK-1 fermion with a KK-1 Higgs boson and a SM fermion 

SM photon with a pair of KK-1 charged Higgs boson 

A pair of KK-1 gauge bosons with a SM or Z° boson 

A pair of KK-1 gluons with a SM gluon 

A pair of SM gluons with a pair of KK-1 gluons 

SM boson with a KK-1 charged Higgs boson and a 

KK-1 pseudoscalar Higgs boson 

SM Z° boson with a pair of KK-1 charged Higgs boson 
SM Z° boson with a KK-1 pseudoscalar Higgs boson and 
a KK-1 scalar Higgs boson 



Table 9: Herwig++ Vertex classes for the UED model. 

HiggsBoundaryMass the value of the Higgs mass at the boundary, the default value is GeV. 
The full list of interfaces for all the classes is provided in the Doxygen documentation. 
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6 Parton Showers 



A major success of the original HERWIG program was its treatment of soft gluon interference 
effects, in particular the phenomenon of colour coherence, via the angular ordering of emissions 
in the parton shower [1,53-61]. Herwig++ simulates parton showers using the coherent branching 
algorithm of [20], which generalizes that used in the original HERWIG program [1-3]. The new 
algorithm retains angular ordering as a central feature and improves on its predecessor in a 
number of ways, the most notable of these being: 

• a covariant formulation of the showering algorithm, which is invariant under boosts along 
the jet axis; 

• the treatment of heavy quark fragmentation through the use of mass-dependent splitting 
functions [62] and kinematics, providing a complete description of the so-called dead-cone 
region. 

In this section we give a full description of the parton shower model and its implementation in 
the program. We begin by introducing the fundamental kinematics and dynamics underlying the 
shower algorithm. This is followed by descriptions of the initial conditions and the Monte Carlo 
algorithms used to generate the showers. Toward the end of the section we discuss how some 
next-to- leading log corrections can be included by a redefinition of the running coupling constant 
and process-specific matrix element corrections. The section concludes with details of the C++ 
code structure. 



6.1 Shower kinematics 

Each colour-charged leg of the hard sub-process is considered to be a shower progenitor. We 
associate a set of basis vectors to each progenitor, in terms of which we can express the momentum 
(<7i) of each particle in the resulting shower as 

q, t = a(p + Pin + q ±i . (6.1) 

This is the well known Sudakov basis. The vector p is equal to the momentum of the shower 
progenitor generated by the prior simulation of the hard scattering process, i.e. p 2 = m 2 , where 
m is the on-shell mass of the progenitor. The reference vector n is a light-like vector that satisfies 
n ■ p > m 2 . In practice n is chosen anticollinear to p in the frame where the shower is generated, 
maximizing n ■ p. Since we almost always generate the shower in the rest frame of the progenitor 
and an object with which it shares a colour line, n is therefore collinear with this colour partner 
object. The q±i vector gives the remaining components of the momentum, transverse to p and n. 
Our basis vectors satisfy the following relations: 

q ±i -p = 0, p 2 = m 2 , q 2 ±i = -q 2 ±i , , g ^ 

q±i • n = 0, n 2 = 0, n-p > m 2 , 

where q_i_j is the spatial component of q±i in the frame where the shower is generated (q 2 ^ > 0). 
Given these definitions, calculating qf, one finds that may be conveniently expressed in terms 
of the mass and transverse momentum of particle i as 

ft = i~f^-^> . (6.3) 
lain ■ p 
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The shower algorithm does not generate the momenta or Sudakov parameters directly. In 
practice what is generated first is a set, each element of which consists of three shower variables, 
which fully parameterize each parton branching. One of these variables parameterizes the scale 
of each branching, the so-called evolution scale, which we shall discuss in more detail below. 
Typically this evolution scale starts at a high value, characteristic of the process that produces the 
progenitors, and continually reduces as the shower develops, via the radiation of particles. When 
the evolution scale has reduced to the point where there is insufficient phase space to produce any 
more branchings, the resulting partons are considered to be on-shell, and the reconstruction of 
the momenta from the shower variables may begin in full. We now define these shower variables. 

The first shower variable we introduce is the light- cone momentum fraction z. Given a 
branching, ij — > i + j n , this parameterizes how the momentum component of the parent parton, 
ij, in the direction of the shower progenitor, is divided between its two daughter partons, i and 
j. We define z as 

z = — = (6.4) 
a rj n-q rj 

For particles in the final state we use a forward evolution algorithm where the parton shower 
consists of a sequence of branchings ij — > i + j, ordered in the evolution scale. For incoming 
particles we use a backward evolution algorithm where we start at the large evolution scale of 
the scattering process and evolve the incoming particles backwards toward the incoming hadron 
to give the mother ij and the sister parton j, again with a decreasing evolution scale. We use 
the definition of z in Eq. (6.4) both for forward and backward parton shower algorithms. 

The second variable used to parameterize a branching is the azimuthal angle, <fr, of the relative 
transverse momentum of each branching p±, measured with respect to the p direction. The 
relative transverse momentum p± is defined to be 

P± = q±i-zq ±rj . (6.5) 

As with the definition of z, this definition of the relative transverse momentum is the same for 
both forward and backward parton-shower evolution algorithms. 

The last, and most important, of the shower variables defining a branching is the evolution 
scale. Parton shower algorithms may be formulated as an evolution in the virtualities of the 
branching partons, or as an evolution in the transverse momentum of the branching products. 
However, a careful treatment of colour coherence effects [1,53-61] reveals that branchings involv- 
ing soft gluons should be ordered in the angle between the branching products. 

The key finding in these studies is that, when soft gluon emissions are considered, branch- 
ings that are not angular ordered do not give any leading logarithmic contributions. This is a 
dynamical effect whereby radiation from the emitting partons, with smaller angular separations, 
interferes destructively in these non-ordered regions. Some intuitive understanding of the effect 
may be gained by considering that a soft gluon, emitted at a large angle from a jet-like con- 
figuration of partons, does not have sufficient transverse resolving power to probe the internal 
jet structure. As a result, it is only sensitive to the coherent sum of the collinear singular con- 
tributions associated with the constituents, resulting in a contribution equivalent to that from 
the original progenitor parton. Destructive interference in the non-ordered region effectively de- 
creases the available phase space for each branching, from the virtuality-ordered region to the 
angular-ordered region. 



L We reserve the tilde notation ij exclusively to denote the parent parton, which decays into daughters i and j. 
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It may be shown that the contributions that angular ordering misses are purely soft and 
suppressed by at least one power of N^, where Nq = 3, the number of colours in QCD. Formally 
then, omitting such contributions amounts to neglecting terms of next-to-leading-log accuracy 
that are also strongly colour suppressed. We stress however, that whereas angular ordering leads 
to an omission of these suppressed higher order terms, other forms of ordering must prove that 
they do not overestimate leading-log contributions. 

For the forward evolution of partons with time-like virtualities, the variable used to achieve 
such ordering, q 2 , is defined according to 



;i-.H 2 = -m! + ^ + -^ (6.6) 

z 1 — z z (1 — z) 



where rrii is the on-shell mass of particle i etc. This definition is arrived at by generalizing 
the FORTRAN HERWIG angular evolution variable, q 2 = q~/ (z (1 — z)), to include the effects 
of the mass of the emitting parton. This may be seen by writing q~ = q-i + qj, and calculating 
<4 {z, p\, q 2 , qf) , which shows 



2 2 

q~ — m~ 
q ' z(l-z) 



(6.7) 



q 2 =m 2 



For showers involving the evolution of partons with space-like virtualities, the evolution variable 
is instead defined by 

ZTTl 2 T) 2 

(1 - z) q 2 = -zm}. + m 2 + —L- - (6.8) 
J 1 — z 1 — z 

Once again this definition of the evolution variable is a generalization of the analogous FORTRAN 
HERWIG angular evolution variable used for initial-state radiation: q 2 = q 2 / (1 — z). Using 

momentum conservation, q^ = q t + qj, we may calculate q 2 q 2 ., q 2 j , whence one finds 



~ 2 = m l ~ Qi 



(6.9) 



9 9 9 9 

g~=m~, q-=mr- 



To see how these variables relate to the angle between the branching products, consider that 
the parton shower is generated in the frame where the light-like basis vector n is anticollinear to 
the progenitor. For forward evolving partons with small time-like virtualities, expanding z and 
g~ in component form, one finds 

2E|(l-cos^)(l + cos^.) 2 
a 2 = — ^ (6 10) 

q (l + COS0i)(l+COS0j) ' 1 J 

where Qi and 6j are the angles between the daughter particles i, j and the progenitor, Qr. is 
the angle between the parent and the progenitor, and 9ij is the angle between the two daugh- 
ters. E-j denotes the energy of the parent. This expression for the time-like evolution vari- 
able in terms of angles is more complicated than the analogous FORTRAN HERWIG formula: 
q 2 = 2Ef. (1 - cos0 y ). This is due to the fact that in FORTRAN HERWIG z was defined to be 
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the energy fraction Ei/En, instead of the light-cone momentum fraction as given in Eq. (6.4). 
Nevertheless, for small angles we find that the Herwig++ and FORTRAN HERWIG evolution 
variables are both given by 

g = ^-(l-0(^)). (6.11) 

When a branching occurs, the daughter partons i and j, with momentum fractions z and 
1 — z, have their starting evolution scales set to zq and (1 — z) q respectively, where zq ~ E$ij 
and (1 — z) q ~ Ej9ij. In this way the maximum opening angle of any subsequent branching is 
8ij, thereby implementing angular ordering. 

For initial-state showers the same QCD coherence argument applies, so in evolving backwards, 
away from the hard process, the angle between the mother of the branching and its final-state 
daughter parton must decrease. Writing the space-like evolution variable (Eq. (6.8)) in terms 
of angles, neglecting parton virtualities, one finds the same form as for the time-like variable 
in Eq. (6.11). This means that once a branching has occurred in the course of the backward 
evolution, the mother of the branching evolves backward from scale q, and the daughter evolves 
forward from scale (1 — z) q, as in the time-like case. 

As stated above, when the evolution in terms of the shower variables has run its course, i.e. 
there is no more phase space available for further emissions, the external particles are taken as 
being on-shell and the reconstruction in terms of the physical momenta can start. First, all of 
the a coefficients in the Sudakov decomposition of each momentum are calculated. This is done 
by first setting a equal to one for final-state progenitors and to the associated PDF light-cone 
momentum fraction x, generated in the preceding simulation of the hard process, for initial-state 
progenitors. Using the defining z relation Eq. (6.4), together with the momentum conservation 
relation = a.{ + Oj, one can iteratively calculate all a values, starting from the hard process 
and working outward to the external legs. 

For final-state showers the q± components of each momentum may be simultaneously calcu- 
lated. Final-state showering cannot change the direction of the progenitor since the transverse 
momentum must be conserved at each branching, hence the q± component of the progenitor 
is zero. The q± components of the branching products are iteratively computed by adding the 
relative transverse momentum, 

P± = (|p_l|cos0, |p_L|sin0,O;O), (6.12) 

to z times the transverse momentum of the mother, to give q±i according to Eq. (6.5); 
q±j = q ± -j — q±i immediately follows by momentum conservation. The magnitude of the relative 

transverse momentum |pj_| = \J—p\ is calculated in terms of the evolution variables z and q 2 
using Eq. (6.6). 

The only remaining Sudakov parameters to be determined are the (3 values. These can be 
obtained once the evolution in terms of the shower variables is complete, by using the fact 
that the external partons are on-shell, in order to compute their 3 coefficients from Eq. (6.3). 
The coefficients of their parent momenta may then be computed using momentum conservation: 
fiu — A + A- The latter step may be iterated until the progenitor is reached, yielding all 3 
coefficients. 

The reconstruction of the initial-state parton showers is slightly different but it follows es- 
sentially the same reasoning. Our aim here has been to simply sketch how the reconstruction 
occurs. More detailed presentations of these procedures will be given later in Sects. 6.4, 6.5 and 
6.6. 
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6.2 Shower dynamics 

With the kinematics defined, we now consider the dynamics governing the parton branchings. 
Each parton branching is approximated by the quasi- collinear limit [62] , in which the transverse 
momentum squared, p^, and the mass squared of the particles involved are small (compared to 



p-ri) but p]_/m is not necessarily small. In this limit the probability of the branching ij 
can be written as 

a s dq 2 



•i.i 



2tt 



i+j 



(6.13) 



where Pfj^j (z, q) are the quasi-collinear splitting functions derived in [62]. In terms of our light- 



cone momentum fraction and (time-like) evolution variable the quasi-collinear splitting functions 
are 
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for QCD and singular SUSY QCD branchings 12 . These splitting functions give a correct physical 
description of the dead-cone region p^ < m, where the collinear singular limit of the matrix 
element is screened by the mass m of the emitting parton. 

The soft limit of the splitting functions is also important. The splitting functions with soft 
singularities P q ^ qg , Pq^qg, Pg^ggi an d Pg^gg, in which the emitted particle j is a gluon, all behave 
as 



lim P~ 



2C r . 



in the soft z 
P 



i _ ml 

q 2 



(6.15) 

1 — 2 

1 limit, where equals Cf for P q ^ qg and Pq^qg, |CV 3 for P g ^ gg , and Ca for 
g -> g g- In using these splitting functions to simulate the emission of a gluon from a time-like 
mother parton ij, associated to a general n parton configuration with matrix element Ai n , one 
is effectively approximating the matrix element for the process with the additional gluon, Ai n+ i, 
by 



+i 



2 2 

q~ — m~ 



\M n \ 



(6.16) 



12 



The Pg^gg splitting presented here is for final-state branching where the outgoing gluons are not identified 
and therefore it lacks a factor of two due to the identical particle symmetry factor. For initial-state branching 
one of the gluons is identified as being space-like and one as time-like and therefore an additional factor of 2 is 
required. 

13 Note that for g — > gg, there is also a soft singularity at z — > with the same strength, so that the total 
emission strength for soft gluons from particles of all types in a given representation is the same: Cf in the 
fundamental representation and Ca in the adjoint. 
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Using the definitions of our shower variables, Eq. (6.4), and making the soft emission approxi- 
mations q^. ~ qi ~ p, qf ~ mf = m~ in Eqs. (6.15, 6.16) we find [23] 

lim 4^T^ r, = -kcasCz (-=- - ^-V. (6.17) 

Recalling that we choose our Sudakov basis vector n to point in the direction of the colour 
partner of the gluon emitter (ij/i), Eq. (6.17) is then just the usual soft eikonal dipole function 
describing soft gluon radiation by a colour dipole [63], at least for the majority of cases where 
the colour partner is massless or nearly massless. In practice, the majority of processes we intend 
to simulate involve massless or light partons, or partons that are light enough that n reproduces 
the colour partner momentum to high accuracy 14 . 

For the case that the underlying process with matrix element M. n is comprised of a single 
colour dipole (as is the case for a number of important processes), the parton shower approxi- 
mation to the matrix element A4 n+ i, Eq. (6.16), then becomes exact in the soft limit as well as, 
and independently of, the collinear limit. This leads to a better description of soft wide angle 
radiation, at least for the first emission, which is of course the widest angle emission in the an- 
gular ordered parton shower. Should the underlying hard process consist of a quark anti-quark 
pair, this exponentiation of the full eikonal current, Eq. (6.17), hidden in the splitting functions, 
combined with a careful treatment of the running coupling (Sect. 6.7), will resum all leading 
and next-to-leading logarithmic corrections [32,64-66]. In the event that there is more than one 
colour dipole in the underlying process, the situation is more complicated due to the ambiguity 
in choosing the colour partner of the gluon, and the presence of non-planar colour topologies. 

In general, the emission probability for the radiation of gluons is infinite in the soft z — > 1 
and collinear q — > limits. Physically these divergences would be canceled by virtual corrections, 
which we do not explicitly calculate but rather include through unitarity. We impose a physical 
cutoff on the gluon and light quark virtualities and call radiation above this limit resolvable. The 
cutoff ensures that the contribution from resolvable radiation is finite. Equally the uncalculated 
virtual corrections ensure that the contribution of the virtual and unresolvable emission below 
the cutoff is also finite. Imposing unitarity, 

V (resolved) + V (unresolved) = 1, (6.18) 

gives the probability of no branching in an infinitesimal increment of the evolution variable dq as 

1 E d ^ -,r ( 6 - 19 ) 

ij 

where the sum runs over all possible branchings of the particle ij. The probability that a parton 
does not branch between two scales is given by the product of the probabilities that it did not 
branch in any of the small increments dq between the two scales. Hence, in the limit dq — > the 
probability of no branching exponentiates, giving the Sudakov form factor 

A(q,q h ) = l[A r ^ tJ (q,q h ) (6.20) 

ij 

14 Even when the colour partner has a large mass, as in e + e~ — > tt, the fact that each shower evolves into the 
forward hemisphere, in the opposite direction to the colour partner, means that the difference between Eq. (6.17) 
and the exact dipole function is rather small in practice. 
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which is the probability of evolving between the scale q^ and q without resolvable emission. The 
no-emission probability for a given type of radiation is 

(S, q h ) = exp {- jT |£ / dz iV.^. . g~') 6 (pi > 0) } . (6.21) 

The allowed phase space for each branching is obtained by requiring that the relative transverse 
momentum is real, or pi > 0. For a general time-like branching ij — > i + j this gives 

z 2 (1 - zf q 2 - (1 - z) m 2 - zm 2 + z {1 - z) m\. > 0, (6.22) 

from Eq. (6.6). 

In practice rather than using the physical masses for the light quarks and gluon we impose a 
cutoff to ensure that the emission probability is finite. We use a cutoff, Q g , for the gluon mass, 
and we take the masses of the other partons to be fi = max (m, Q g ), i.e. Q g is the lowest mass 
allowed for any particle. 

There are two important special cases. 

1. q — > qg, the radiation of a gluon from a quark, or indeed any massive particle. In this case 
Eq. (6.22) simplifies to 

z 2 (l - zfq 2 > (1 - zf /i 2 + zQ 2 g , (6.23) 
which gives a complicated boundary in the (q, z) plane. However as 

(1 _ z f ^ + ^ > (1 _ z) 2 ^ z 2 Q 2 (6 24) 

the phase space lies inside the region 

^ < z < 1 - ^ (6.25) 

q q 

and approaches these limits for large values of q. In this case the relative transverse 
momentum of the branching can be determined from the evolution scale as 



p± = V(l - z) 1 (z^ - v?) - zQl (6.26) 

2. g — > gg and g — > gg, or the branching of a gluon into any pair of particles with the same 
mass. In this case the limits on z are 



z_ <z < z+, z± = ^\ l±^l - and q > 4/i. (6.27) 

Therefore analogously to Eq. (6.25) the phase space lies within the range 

^ < z < 1 - ^. (6.28) 

q q 

In this case the relative transverse momentum of the branching can be determined from 
the evolution scale as 

p ± = Jz 2 (l-zfq 2 -fi 2 . (6.29) 
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These two special cases are sufficient for all the branchings currently included in the simulation, 
although the general case of three unequal masses for the particles in the branching is supported. 

The cutoff parameter, Q g , is the minimum virtuality of the gluon. However, if we consider 
the phase space that is available to the parton shower we would expect a natural threshold of 
order m + Q g for gluon emission from a quark of mass m. In practice for the radiation of a gluon 
from a quark, Eq. (6.26) gives a threshold that behaves as Qthr — 1-15 (m q + 2Q g ). This means 
that the phase-space limit is well above our expectation, particularly for heavy quarks. 

There is no reason why Q g should be the same for all quark flavours. Therefore, we have 
chosen to parameterize the threshold for different flavours as 



where a [aParameter=0.3] and b [bParameter=2.3] are parameters chosen to give a threshold 
Qthr — P m q + 8, with (3 = 0.85, in order to slightly reduce the threshold distance for heavier 
quarks. As a result, the threshold for radiation from heavy quarks is closer to its physical limit. 
The parameter 5 is tuned to data as [cutoffKinScale=2.8 GeV] and, only relevant for partons 
heavier than the bottom quark, the parameter c is chosen to prevent the cutoff becoming too 
small, [cParameter=0.3 GeV]. 

The formalism discussed above allows us, if given a starting scale %, to evolve a parton 
down in scale and generate the next branching of this particle at a lower scale. The no-emission 
probability encoded in the Sudakov form factor is used to generate (q, z) for this branching. This 
procedure can then be iterated to generate subsequent branchings of the particles produced until 
no further emission occurs above the cutoff. 

6.3 Initial conditions 

Before we can simulate possible radiation from a hard process we need to know the initial con- 
ditions, i.e. the scale qh from which to start the evolution. The initial conditions for the parton 
shower are determined by the colour flow in the hard process [3]. For each particle involved in the 
hard process a colour partner is chosen. In the case of particles in the fundamental representation 
of the SU (3) gauge group this choice is unique, at least in processes where baryon number is 
conserved. In the case of a gluon a uniform random choice is made between the two possible 
partners. In processes involving baryon number violation a uniform random choice is made be- 
tween all the potential colour partners [67,68]. The direction of this colour partner determines 
the maximum angle for emission of QCD radiation from a particle in the angular-ordered parton 
shower. 

Following the choice of the colour partner the maximum scale for radiation from the particle 
must be calculated, as must the choice of the p and n reference vectors defined in Eq. (6.1). We 
always take the choice of p along the direction of the radiating particle but the choice of n is 
related to the direction of the colour partner. 




(6.30) 
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6.3.1 Final-final colour connection 

The easiest case to consider is the colour connection between two final-state particles, b and c. 
Working in their centre-of-mass frame, we may write their momenta as 

Pb = l -Q (0, A; 1 + b - c) Pc = X -Q (0, -A; 1 - b + c) , (6.31) 

where Q 2 = (p b +p c ) 2 , b = ml/Q 2 , c = m 2 c /Q 2 and 

A = A (1, b, c) = Vl + b 2 + c 2 - 2b - 2c - 2bc (6.32) 

is the Callan function. 

In order that the soft region of phase space is fully covered, the initial evolution scales for b 
and c {qubi Qhc) are related by 

(k b -b)(k c -c) = ^(l-b-c + \) 2 , (6.33) 

where k b = q 2 lb /Q 2 , k c = ql c /Q 2 [20]. By varying the starting scales of the individual particles 
we can control how much radiation is generated from each of them, in order to assess the uncer- 
tainties. In practice we currently allow four choices controlled by the FinalFinalConditions 
switch: 

Symmetric The most symmetric choice of the initial conditions, giving equal amounts of 
radiation from both partons is given by 

k b = I(i + 6- c + A), k c = I(l-6 + c + A). (6.34) 

This is our default choice [FinalFinalConditions= Symmetric] . 

Coloured The largest emission scale that is possible for radiation from one of the particles is 
given by 

K b = 4(l-2Vb-b + c\ . (6.35) 

The [FinalFinalConditions= Coloured] choice of initial conditions maximizes the initial evo- 
lution scale for the shower of the coloured particle. Naturally, this therefore minimizes the phase 
space volume available for the first emission from the anti-coloured parton. 

AntiColoured The [FinalFinalConditions=AntiColoured] choice of initial conditions is 
the converse of the [FinalFinalConditions=Coloured] choice. 

Random Selecting the option [FinalFinalConditions=Random] , the program randomly 
sets the initial evolution scales according to the Coloured or AntiColoured options, for each 
final-state pair of colour partners, for each event. 

As stated in Sect. 6.1 the p basis vector (Eq. (6.1)) is given by the momentum of the progenitor 
as it was generated in the initial simulation of the hard process. The light-like basis vector n is 
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chosen to be collinear with the colour partner in the rest frame of the coloured connected pair, 
i.e. in simulating radiation from b, n is defined to be 

n — -Q (0, —A; A) . (6.36) 

To simulate parton showering from c, we simply reverse the spatial components of n in Eq. (6.36). 

6.3.2 Initial-initial colour connection 

Here again we opt to work in the rest frame of the colour partners, so that the momenta of the 
particles are 

p 6 = ig(0,l;l) p c =i<2(0,-l;l), (6.37) 

where Q is the partonic centre-of-mass energy of the collision. 

In this case the requirement that the soft region of phase space is smoothly covered is simply 

kbk c = 1. (6.38) 

Contrary to the case of the final-final colour connection, there is no upper bound on the values 
of Rf, or k c , i.e. there is no choice that maximizes the phase space available to one parton relative 
to the other (at least none that might reasonably be expected to give sensible results). Currently 
only the most symmetric choice is implemented, i.e. K\ ) = k c = 1. 

In this case, as we assume that the incoming particles are massless, we can simply take the p 
reference vector to be the momentum of the beam particle from which the emitting parton was 
extracted and the n reference vector to be the momentum of the beam particle from which its 
colour partner was extracted. The fact that p is parallel to the momentum of the emitting parton 
makes it easier to reconstruct the momenta of the shower particles in terms of the fraction of the 
beam momentum they carry. 

Finally, defining the p and n vectors as being equal to the beam momenta rather than the 
actual parton momenta does not affect our earlier assertions relating to the soft limit of the 
splitting functions, since Eq. (6.17) is clearly invariant under overall rescalings of the dipole 
momenta n and p. 

6.3.3 Initial-final colour connection in the hard process 

Consider the initial-final-state colour connection in the context of a process a + b — > c, where a 
is a colour-singlet system and b and c are colour connected, e.g. deep inelastic scattering. As in 
the last two cases we work in the rest frame of the colour dipole, in this case the Breit frame, 
where we may write 

p b = ~Q(0,l + c;l + c), p c = Iq (0,-1 + c;l + c), (6.39) 

with Q 2 = —pi. 

To achieve a smooth matching of the phase space for the first emission from parton 6's shower 
with that of parton c's shower, at wide angles, requires the initial evolution scales (qhb, Qhc) to 
obey 

k h {k c -c) = {l + cf . (6.40) 
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In practice, we opt to assign more-or-less the same phase space volume to each shower, i.e. we use 
the most symmetric choice: k b = 1 + c, k c = 1 + 2c. Of course, a larger or smaller combination 
that satisfies Eq. (6.40) is also allowed. 

For emission from the final-state particle, the p vector is taken to be the momentum of 
the radiating particle and the n reference vector is set equal to the momentum of the beam 
particle from which the initial-state colour partner was extracted. For emission from the initial- 
state particle the p vector is defined to be the momentum of the beam particle from which the 
radiating parton was extracted and 

n=~Q(0,-l-c;l + c), (6.41) 

in the Breit frame. As discussed at the end of the description of the initial-initial colour con- 
nection, the normalization of n and/or p, does not affect the eikonal dipole limit of the splitting 
functions Eq. (6.17). 

6.3.4 Initial-final colour connection in decays 

The Herwig++ shower differs from other approaches in including initial-state radiation from a 
decaying coloured particle, as well as final-state radiation from the coloured decay products. This 
is required in order to ensure that the full soft region of phase space is filled by radiation from 
the parton shower [20,23]. 

Consider the decay b — > ac, where b and c are colour partners and a is a colour singlet system, 
in the rest frame of the decaying particle. In this frame the momentum of b and its colour partner 
c are 

p b = m b (0, 0; 1) , p c = ^m b (0, A; 1 - a + c) , (6.42) 

where c = m^/m^ and hence A = A(l, a, c) where a = m^/m^. 

In this case the requirement that the full soft region of phase space is filled by radiation from 
the parton shower gives 

(«6-1)(k c -c) = i(l-a + c + A) 2 . (6.43) 

While there is no limit on the value of k b as with the final-final colour connection the maximum 
value of k c is 

k c = 4 (1 + a - 2^fc - c) . (6.44) 

We support three choices for the values of the scales controlled by the switch 
InitialFinalDecayConditions 

Symmetric The most symmetric choice of initial conditions is 

k b = ^(3-a + c + A), k c = - (1 - a + 3c + A) , (6.45) 

which is the default choice [InitialFinalDecayConditions=Symmetric] . 
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Maximal The maximal choice corresponds to generating the maximal amount of radiation 
from the final-state particle, i.e. k c is given by Eq. (6.44). This corresponds to 
[InitialFinalDecayConditions=Maximal] . 



Smooth In this case the initial conditions are chosen in order to guarantee that, in addition 
to covering the full soft region, the radiation pattern smoothly changes between the region filled 
by radiation from b and c. In this case 

2A 

Kb = \ — T\ r^~> (6 - 46) 

A - (1 - y/c) + a 

with k c obtained from Eq. (6.43). This option is obtained by setting 

[InitialFinalDecayConditions=Smooth] . In, for example, top decays, this choice leads to 
more radiation from the decaying particle and less from its colour partner than either of the 
other options 15 . 

For radiation from the decaying particle, p is chosen to be the momentum of the decaying 
particle and 

n=im 6 (0,l;l), (6.47) 

in its rest frame, i.e. n is aligned with the colour partner. 

In the case of radiation from the final-state particle, p is set equal to its momentum, as 
generated in the hard decay process, however, there is no obvious choice of n related to the 
colour partner, since we are working in its rest frame. We therefore choose n such that it is in 
the opposite direction to the radiating particle in this frame, i.e. 

n = i(0,-A;A). (6.48) 

A more rigorous approach to this problem was carried out in [23], using a more generalized 
splitting function, derived assuming a massive gauge vector n. This feature is not implemented 
in the standard released code, since any related deficiency in the shower is completely avoided 
by using the associated matrix element correction (Sect. 6.8). 



6.4 Final-state radiation 
6.4.1 Evolution 

The parton shower algorithm generates the radiation from each progenitor independently, modulo 
the prior determination of the initial evolution scale and the n and p basis vectors. Consider 
then, the evolution of a given final-state progenitor, downward from its initial evolution scale 
Given that A (q, qh) gives the probability that this parton evolves from scale q\ to q without any 
resolvable branchings, we may generate the scale of this first branching (q) by solving 

A(q,q h )=TZ, (6.49) 

15 In the extreme limit c — > 0, e.g. if in top decays the bottom quark is considered massless relative to the top, 
Kb — > oo and k c — > 0, meaning that emission only comes from the decaying top quark and none at all from the 
massless bottom quark. This is because in the limit of a massless bottom quark radiation from the top quark 
gives the correct dipole distribution in the soft limit. 
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where TZ is a random number uniformly distributed between and 1. 

In the FORTRAN HERWIG program this equation was solved by a brute force numerical 
calculation, using an interpolation table for A (g, In Herwig++ an alternative approach is 
used, which determines the scale of the branchings without the need for any explicit integration 
of the Sudakov form factor [69] . The method involves generating each branching according to a 
crude Sudakov form factor, based on an overestimated branching probability (Eq. (6.13)), simple 
enough that Eq. (6.49) can be solved analytically. Each of these crudely determined branchings 
is subject to a vetoing procedure based on a series of weights relating to the true form factor. 
In this way the overestimated, crude emission rate and emission distribution is reduced to the 
exact distribution. 

The first ingredient we need in order to implement the algorithm is therefore a crude approx- 
imation to the Sudakov form factor (Eq. (6.21)), for each type of branching of a parent parton 
ij, ij ~ > i + j- We write these as 



Qh 

over, res. 



where 

do 2 r z °+ er ry OVCT 

dp ovcr,res. = ^_ dz pover / n (g^) 

is the overestimated probability that a resolvable branching occurs in the interval [q 2 , q 2 + dq 2 }. 
Overestimates of the splitting functions and the coupling constant are denoted P? Y ^ (z) > 
Pfj^ij (z,q) and a| ver > as(z,q), while the limits z°± er also denote overestimates of the true 
z integration region 16 for all q. To solve Eq. (6.49) analytically we also require that Pf™^ {z) 
should be analytically integrable and, in order to generate z values, this integral should be an 
invertible function of z. 

Using this simplified Sudakov form factor we may analytically solve 

q 2 = q 2 h K$, (6.52) 

where 

^pOver,res. 

which can be thought of as the number of emissions per unit of the shower formation 'time' lng 2 
(for the crude distribution this is a constant). It is clear from Eq. (6.52) how increasing this rate 
r causes the first branching to be generated 'sooner', closer to q^. When a value is generated for 
the evolution scale of a branching, an associated z value is then generated according to 

z = r 1 [I (z°_ ver ) + n' (I (z° + ver ) - I (z°_ ver ))} , (6.54) 

where / (z) is the primitive integral of PfJ ev i; j i z ) over z, I" 1 is the inverse of / and TZ' is a 
uniformly distributed random number in the interval [0, 1]. 
We now reject these values of q and z if: 



16 



The overestimates of these limits were given in Eqs. (6.25,6.28). 
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the value of z lies outside the true phase-space limits, i.e. if p 2 ± < 0; 



as(z,q) -i-, . 

a s 

w pover ( z \ ~- /v 2) 

where are random numbers uniformly distributed between and 1. 

If we reject the value of q we repeat the whole procedure with qt = q until either we accept a 
value of q, or the value drops below the minimum value allowed due to the phase-space cutoffs, 
in which case there is no radiation from the particle. As shown in [69] this procedure, called the 
veto algorithm, exponentiates the rejection factors and generates the values of q and z according 
to Eq. (6.49) for one type of branching. 

This procedure is repeated to give a value of the evolution scale for each possible type of 
branching, and the branching with the largest value of q is selected, which then generates both 
the type of branching, its scale, and the momentum fraction according to Eq. (6.49), as required. 

The relative transverse momentum for the branching p± (Eq. (6.5)) is then calculated, using 
Eq. (6.26) or Eq. (6.29) depending on the type of branching. Currently the azimuthal angle of p± 
is randomly generated between and 2ir about the direction of the progenitor (the Sudakov basis 
vector p), although in future this will change when we include spin correlations in the parton 
shower as described in [27-30]. 

The requirement of angular ordering, as discussed in Sect. 6.1, defines the initial scales for 
the daughter particles, %i and qhj, produced in each branching, ij — > i + j, to be 

q~hi = zq, q hj = (1 - z) q, (6.55) 

where q and z, are the evolution scale and light-cone momentum fraction of the branching. By 
imposing these upper bounds on the evolution scale of the emitted partons, subsequent branchings 
will have a nesting of the angular separation of the resulting daughters, where each one is smaller 
than the one preceding it. 

All of the steps above, required to generate the shower variables associated with this initial 
branching, may then be repeated for the daughter partons, and their daughter partons, should 
they also branch. All showering terminates when the evolution scale (q) for each final-state parton 
falls below its minimum value, when there is no phase space for any more resolvable emissions. 
The resulting partons, at the end of each shower, are deemed to be on constituent mass-shell, as 
defined in Sect. 7, at which point the perturbative parton shower evolution is no longer sensible, 
since hadronization effects dominate at these scales. 

6.4.2 Kinematic reconstruction 

At this point we have a set of partons produced in the parton shower from each of the progenitor 
partons, the scales q at which they are produced, the momentum fractions z and azimuthal angles 
cf) of the branchings. Mapping these kinematic variables into physical momenta is what we call 
kinematic reconstruction. We will now describe this procedure for showers generated by final- 
state progenitors. First, the kinematics of the individual showers are reconstructed by putting 
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the external masses on their constituent mass-shell 17 and working back through the shower, as 
described in Sect. 6.1. 

The shower evolution causes all progenitor partons, J, produced in the hard process to gain 
a virtual mass, i.e. the progenitor partons, which initiated the jets, are no longer on mass shell, 
q 2 7^ m 2 . We want to preserve the total energy of the system in the centre-of-mass frame of 
the hard collision. If the momenta of the progenitor partons before the shower evolution were 

Pj = ^pj; -s/pj + in this frame, then 

n 

j=i 

while the sum of the spatial momenta is zero. As the jet parents have momenta 

qj = ^qj; -\/q j + qjj after the parton showering, we need to restore momentum conservation 

in a way that disturbs the internal structure of the jet as little as possible. The easiest way to 
achieve this is by boosting each jet along its axis so that their momenta are rescaled by a common 
factor k determined from 

n 

j=i 

which can be solved analytically for two jets or numerically for higher multiplicities. For every 
jet a Lorentz boost is applied such that 

QJ = (<Lrl yj<& + <#) ^ QJ = (kpr, y/k^j + . (6.58) 

Applying these boosts to each of the jets, in the centre-of-mass frame of the collision, ensures 
global energy-momentum conservation. Typically the rescaling parameters k are close to unity, 
hence the resulting boosts and rotations are small. 



6.5 Initial-state radiation 
6.5.1 Evolution 

As stated in Sect. 6.1, in generating the initial-state radiation we use a backward evolution 
algorithm, starting with the space-like daughter parton that initiates the hard scattering process, 
i, and evolving it backward to give its space-like parent, ij, and time-like sister parton j. This 
evolution algorithm therefore proceeds from the high scale of the hard process to the low scale 
of the external hadrons. Such a procedure is greatly more efficient than the alternative forward 
evolution algorithm, which would start from the incoming beam partons and evolve them to the 
scale of the hard collision. This is because the forward evolution cannot be constrained to end on 
the x and Q 2 values associated to the hard process, which in turn makes it impossible to perform 
importance sampling of any significant resonant contributions. 

While forward evolution would dynamically generate the parton distribution functions (PDFs), 
backward evolution uses the measured PDFs to guide the evolution. As with the final-state 

17 Thc Herwig++ shower allows these masses to be set to zero so that an alternative hadronization model, rather 
than the cluster model, can be used. 
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shower, the initial conditions for the initial-state shower are determined by the colour partners 
of the incoming particles (Sect. 6.3.2). 

The angular-evolution variable q 2 for space-like showers was defined in Eq. (6.8). We shall 
work exclusively with light initial-state partons so we take mp — rrii — 0, and rrij = fi if j is a 
quark and rrij = Q g if j is a gluon, to regulate the infrared divergent regions, hence Eq. (6.8) 
simplifies to 



„ zm 2 - p 2 ± 

^^tf' (6 - 59) 



where p]_ = —p 2 ± (Eqs. (6.5,6.12)). 

From the requirement that > 0, Eq. (6.59) implies an upper limit on z, 



2g 2 V V 2g 2 



z<z+ = l + ^-\ [1 + ^] -I. (6.60) 



In addition, if the light-cone momentum fraction of parton % is x, we must have z > x to prevent 
the initial-state branching simulation evolving backward into a parent with x > 1. 
In this case the Sudakov form factor for backward evolution is [3,70] 

A (x, q, q h ) = Y[ ^ij^ij ( x ' Qf Qh) » ( 6 - 61 ) 

where the Sudakov form factor for the backward evolution of a given parton i is 

A 5 - y (*, J, & ) = exp { - jT f f d, SO §gje (pi > 0) } (6.62) 

and the product runs over all possible branchings ij — > i+j capable of producing a parton of type 
%. This is similar to the form factor used for final-state radiation, Eq. (6.21), with the addition 
of the PDF factor, which guides the backward evolution. 

The backward evolution can be performed using the veto algorithm in the same way as the 
forward evolution. We need to solve 

A(x,q,q h )=K, (6.63) 

to give the scale of the branching. We start by considering the backward evolution of % via a 
particular type of branching, ij — > i + j. We can obtain an overestimate of the integrand in the 
Sudakov form factor 

f r« h c\fi' 2 r z + er ry OVQT } 
Ag- . q, ~q h ) = exp | - J J dz -j- P^.. (z) PDF ovcr (z) j , (6.64) 

where Pp™.. (z) , a^ ei and the overestimate of the limits must have the same properties as for 
final-state branching. In addition 

PDF ovcr (z) > ~ zf l j r^ \/z, q, x. (6.65) 
xfi [x, q) 
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In this case the product -P-™^- ( z ) PDF ovcr (z) must be integrable and the integral invertible. If 
we define 



„over — over 



r = ^trj dz (*) PDF ° vcr (*) > ( 6 - 66 ) 

then we can solve Eq. (6.62) using this overestimated Sudakov giving 

q 2 = q 2 h rt. (6.67) 

The value of z can then be generated according to 

z = r 1 [I ( x ) + n 1 (I (z° + ver ) - I (z))] , (6.68) 

where I (z) = J dz P~ VCT . . (z)PDF over (z), J -1 is the inverse of / and 1Z' is a random number 
uniformly distributed between and 1. 
We now reject these values of q and z if: 

• the value of z lies outside the true phase-space limits, i.e. if < 0; 



-..over 



P- ..(z,q) ^ 
povcr ( z \ ■■ '^2; 
Ij— »ij 



< Ha 



PDF over (a) 

where 72-1^,3 are random numbers uniformly distributed between and 1. 

As with the final-state branching algorithm, if a set of values of q and z, generated according to 
the approximate form factor in Eq. (6.64) is rejected, a further set is then generated by repeating 
the process with q h = q in Eq. (6.64). This procedure is repeated until either a generated set of 
branching variables passes all four vetoes, or the generated value of q falls below the minimum 
allowed value, in which case the showering of the particle in question ceases. To determine the 
species of partons involved, a trial value of q is generated for each possible type of branching 
and the largest selected. By applying the four vetoing criteria to each emission generated by the 
approximate, overestimated, Sudakov form factor, the accepted values of q and z are distributed 
according to the exact Sudakov form factor, Eq. (6.62) [69]. 

When a branching is generated, the relative transverse momentum p± (Eqs. (6.5, 6.12)) is 
calculated according to Eq. (6.59). At present the azimuthal angle associated to each p± is 
randomly generated between and 2ir, although in future this will change when we include spin 
correlations in the parton shower as described in [27-30]. In the case of backward evolution the 
angular ordering requirement is satisfied by simply continuing the backward evolution downward 
in q, starting from the value generated in the previous generated branching. 

As stated above, when the evolution scale has reduced to the point where there is no more 
phase space for further resolvable branchings, the backward evolution ends. The incoming particle 
produced in the last backward branching, assumed to be on-shell (massless), has no transverse 
momentum, since this is measured with respect to the beam axis 18 . This final parton also has a 



18 HerwigH — h supports the option of including a non-perturbative intrinsic transverse momentum for the partons 
inside the incoming hadron, as described in Appendices B.9 and D, which can give the initial incoming parton a 
transverse momentum. 
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light-cone momentum fraction x/ Yli z ii with respect to the incoming hadron's momentum, where 
x is the light-cone momentum fraction generated in the initial simulation of the hard process, 
and the product is comprised of all z values generated in the backward evolution. 

Before any momentum reconstruction can begin, we must simulate the effects of final-state 
showers from each time-like daughter parton j, generated from the backward evolution of each 
space-like parton i, in branchings ij — > i + j. As discussed in Sect. 6.1, for such a branching 
occurring at scale q with light-cone momentum fraction z, angular ordering is achieved by evolving 
j down from an initial scale q^ — (1 — z) q. This initial condition ensures that for each parton j, 
the angular separation of any of fs subsequent branching products is less than the angle between 
j and j's sister %. 

This algorithm is all that is needed to generate the values of the scales, momentum fractions 
and azimuthal angles, for the evolution of both the incoming particles and the time-like particles 
emitted in their backward evolution. These values are sufficient for us to determine the momenta 
of all of the particles in the associated showers, to perform the kinematic reconstruction. 

6.5.2 Kinematic reconstruction 

The kinematic reconstruction begins by finding the last initial-state particle produced in the 
backward evolution of each of the beam particles. This parton's momentum is calculated as 
described in the previous section. The momentum of the final-state time-like jet that it radiates 
is then reconstructed in the same way as for the final-state shower. Knowing the momenta of 
the former light-like parent parton and the latter final-state, time-like, daughter parton, the 
momentum of the initial-state, space-like, daughter, follows by momentum conservation. This 
process is iterated for each initial-state branching, eventually giving the momentum of the space- 
like progenitor parton, colliding in the hard process. 

The reconstructed momentum of the colliding parton incident from the +z direction is denoted 
q & , and that of the colliding parton incident from the —z direction is denoted q e . 

The final reshuffling of the momentum then proceeds in different ways depending on whether 
the colour partner is an initial- or final-state parton. 

Initial-State partner As discussed in Sect. 6.3.2 the hadronic beam momenta, p© and p e , 
then define the Sudakov basis for the initial-state shower algorithms, in terms of which we have 



The Sudakov coefficients may be calculated using the fact that and p e are light-like and 
orthogonal to the qj_ component: 



where s = 2p e ■ p e , the hadronic centre-of-mass energy squared. The q± components follow by 
subtracting + P®P® from the reconstructed momentum q^. 

Through the emission of initial-state radiation the colliding partons acquire both space-like 
virtualities and transverse momenta, of which they had neither in the initial simulation of the 
hard process. Consequently, whereas momentum conservation in the prior simulation of the hard 
process implies that the total initial- and final- state momentum are equal to p cms = x e p e +x e p e , 
we now have a momentum imbalance between the two: g e + q e ^ x e p e + x Q p e . 



q® = a@p© + P®p® + q±®- 



(6.69) 



«© = 2 p® ■ q®/ s 



P® = 2p@ • q®/s, 



(6.70) 
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In order to return to a momentum conserving state we choose to rescale the energies and 
longitudinal momenta of the colliding initial-state partons, in a way that preserves the invari- 
ant mass and rapidity of the centre-of-mass system. The transverse momentum of the emitted 
radiation can only be absorbed by the final-state system. When the rescaling factors have been 
determined, we can then calculate a Lorentz boost that produces the same effect. This boost 
can then be applied to all elements of the initial-state shower, including the final-state jets they 
emit. 

The energies and longitudinal momenta of the colliding partons are rescaled by two factors, 
fc© and k e , giving shuffled momenta q'^ and q' e , according to 

<?@ = a@k@p@ + j^p® + q±®- (6.71) 

In simulating the hard process the momentum of the partonic centre-of-mass system was given 
by 

Pcms = 2;ePe + x e p e (6.72) 
and in terms of the shuffled momenta it is 

<?cms = (^®k e + p ffi + (a e k e + y^j p e + q ±(B + q ±e . (6.73) 

Imposing that the centre-of-mass energy generated in the simulation of the hard process is 
preserved, q^ s = pl ms , the Sudakov decompositions of Eqs. (6.72,6.73), imply that the rescalings 
fc® and k e obey the relation 

a®a Q s kl Q + ((a©/?© + a e (3 e - x e x e ) s + (q ±& + q± e f) k &e + /? e /?eS = 0, (6.74) 

where fc© e = k^k e . The further imposition that the rapidity of the partonic centre-of-mass is 
preserved requires that the ratio of the p e coefficient to the p e Sudakov coefficient in q' cms should 
equal that in p cms . This implies a second constraint on fc® and k e 

, 2 _ , x e + «e^®e (f{ 7 c\ 

fc ffi- fc ©e — — 7 -~7T- y°-' b > 

Xq « e K, ffiQ -f- Pq 

The two relations in Eqs. (6.72,6.73) fully determine the k^ and k e rescaling factors. Having 
solved these equations for and k e we go on to determine a longitudinal boost for each initial- 
state jet such that 

q @ q' @ . (6.76) 

This boost may then be applied to all elements of the initial-state shower including any final-state 
partons/jets that they emit. 

This procedure is sufficient for the production of colour-singlet systems, such as electroweak 
gauge bosons in the Drell-Yan process. 

Final-State partner For systems that have an initial-state parton that is colour connected to 
a final-state parton the reconstruction is performed in their Breit frame in order to preserve the 
Q 2 of the system in, for example, DIS processes. 
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The momenta of the initial- and final-state jets are first reconstructed as described above for 
initial-state jets and in Sect. 6.4.2 for final-state jets. The momenta of the jet progenitors which 
are now off-shell are then boosted to the Breit-frame of the original system before the radiation. 
We take pb to be the momentum of the original incoming parton and p c to be the momentum of 
the original outgoing parton and p a = p c — pb, therefore in the Breit-frame 

Pa = Q(0, 0,-1,1). (6-77) 

We can then construct a set of basis vectors, similar to the Sudakov basis defined in Sect. 6.3.3 
for the initial-final colour connection, 

m = Q(0,0,l;l), n 2 = Q(0, 0, -1; 1). (6.78) 

The momenta of the off-shell incoming parton can then be decomposed as 

Qm = ttin^l + An™2 + Q±, (6.79) 

where a in = ^^ L , (3 = and q± = q m — a m rii — An^2- In order to reconstruct the final-state 
momentum we first apply a rotation so that the momentum of the outgoing jet is 

g ou t = «out^i + A>ut™2 + q±, (6.80) 

where A>ut is taken to be one and the requirement that the virtual mass is preserved gives 
cw = 9 ° ut+p± where q\ = —p\. The momenta of the jets are rescaled such that 



<?in,out = «in,out^in,out^l + J^ L U 2 + <?_L, (6.81) 



Wn,out 

which ensures the virtual mass of the partons is preserved. The requirement that the momentum 
of the system is conserved, i.e. 

Pa = g out -gL = Q(0,0,-l;0), (6.82) 

gives 

c^m^m ~ QWt&out = ~ ■ (6.83a) 

An Pout 1 fn 001 \ 

k~~V~ = "2- ( } 

•"in "'out ^ 

As with the initial-initial case once the rescalings have been determined the jets are transformed 
using a boost such that 

<?in,out ?L,out- (6-84) 

The procedures described above are sufficient for simple cases such as the Drell-Yan produc- 
tion of vector bosons in hadron-hadron collisions or deep inelastic scattering. In general however 
the colour structure of the event, particularly in hadron collisions requires a more general proce- 
dure. 

In general, from Herwig++ version 2.3, the following procedure is used to reconstruct the 
kinematics of the hard process. First the colour structure of the hard process is use to construct 
colour singlet systems from the jet progenitors. Depending on the result different approaches are 
used. 
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• If the incoming particles are colour neutral then any final-state colour singlet systems are 
reconstructed as described in Sect. 6.4.2, for example in e + e~ —> qq. 

• If there is a colour- singlet system consisting of the incoming particles together with a 
number of final-state colour singlet systems, e.g. Drell-Yan vector boson production, then 
the kinematics are reconstructed as described above for the initial-initial system. The 
final-state systems are then reconstructed in their rest frames as described in Sect. 6.4.2 
and boosts applied to ensure the recoil from the initial-state radiation is absorbed by the 
final-state systems. 

• If the system consists of colour-neutral particles and an initial-final state colour connected 
system, e.g. deep inelastic scattering, then the kinematics are reconstructed as described 
above for an initial-final system. 

• If the system consists of two separate initial-final state colour connected systems together 
with a number of colour-singlet final-state systems, for example Higgs boson production via 
vector boson fusion or qq — > ti, then the colour- singlet initial-final systems are reconstructed 
as described above and the final-state systems as described in Sect. 6.4.2. 

• In general in hadron-collisions the hard process cannot be decomposed into colour singlet 
systems and a general procedure which preserves the rapidity and mass of the hard collision 
is used. The initial-state jets are reconstructed as discussed above for the initial-initial 
connection. The final-state jets are then reconstructed in the partonic centre-of-mass frame 
of the original hard scattering process as described in Sect. 6.4.2. This is effectively the 
same as reconstructing them in the q' cms rest frame, since the kinematic reconstruction for 
initial-initial connection preserves the invariant mass of the hard process. In the end, the 
jets originating from the final-state particles in the hard process are boosted back to the 
lab frame, where they have a total momentum q' crns . 

This procedure uses the underlying colour flow in the hard process to determine how global energy 
and momentum conservation is enforced where possible and resorts to the general approach used 
before Herwig++ 2.3 when this is not possible. It is still possible to use the general procedure 
which ignores the colour flow for all process using [ReconstructionOption=General] switch 
rather than the default option which uses the colour structure where possible 
[ReconstructionOption=Colour] . 

6.5.3 Forced splitting 

After the perturbative shower evolution has terminated, the cluster hadronization model may 
necessitate some additional forced splitting of the initial-state parton that results. In hadronic 
collisions we require the external initial-state partons, which give rise to the first hard interaction, 
to be valence quarks (antiquarks), colour triplet states. This allows us to treat each proton 
(antiproton) remnant as a diquark (antidiquark) which will be in a colour antitriplet/triplet 
state, in order to keep the incoming hadron colour neutral. Modelling the dissociation in this 
way allows for a simple, minimal, hadronization of the remnant in the cluster hadronization 
model. 

Usually, the perturbative evolution, which is guided by the PDFs, will terminate on a valence 
quark, since the evolution works towards large x and small Q 2 . In the cases where this has not 
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happened, we force the resulting initial-state parton to undergo one or two additional splittings. 
The generation of these additional forced splittings is largely based on the same principles as 
that of the perturbative splittings. 

In the perturbative evolution the scale of the PDFs is frozen at a value Q s for values Q < Q s . 
The default value of Q s is chosen to be small, close to the non-perturbative region but still above 
typical values for the parton shower cutoff ([PDFFreezingScale=2.5*GeV]). This freezing 
scale leaves a little phase space for the (non-perturbative) forced splittings. The forced splittings 
are generated in much the same vein as the perturbative splittings. The evolution starts at Q s 
and the next branching scale is distributed according to dQ/Q, with a lower limit determined by 
the available phase space. The z values are determined from the splitting functions in the same 
way as in the perturbative evolution. The splittings are reweighted by ratios of PDFs as in the 
perturbative evolution. There is only one slight difference, the evolution of the PDFs themselves 
with Q is frozen below Q s . Nevertheless, this reweighting gives the right flavour content of the 
initial hadron. For example in the case of a proton we produce twice as many u quarks as d 
quarks. To force the evolution to end up on a valence quark, we only allow one or two flavours 
in the evolution: 

1. If the initial parton is a seaquark (q) or -antiquark (g), it is forced to evolve into a gluon, 
emitting a q or q, respectively. 

2. If the initial parton is a gluon, from either the perturbative evolution or the forced splitting 
of a seaquark, it is forced to evolve into a valence quark, emitting the matching antiquark. 

In the initial-state showering of additional hard scatters we force the backward evolution of 
the colliding partons to terminate on a gluon. We therefore only need the first kind of forced 
splitting in this case. This gluon is assumed to be relatively soft and branches off from the 
remnant diquark. Again, this allows us to uniquely match up the final-state partons to the cluster 
hadronization model. We should note that the emitted partons from these forced splittings, as 
well as the remnant diquarks, will show up in the event record as decay products of a fictitious 
remnant particle, in order to distinguish them from those which originate from the perturbative 
evolution. Additional details about the colour structure and the event record can be found in [8] . 

6.6 Radiation in particle decays 

In general the hard processes simulated by Herwig++ consist of 2 — > n scatterings. These are 
generated by first using the relevant matrix elements to produce an initial configuration, and 
then initiating parton showers from the external legs. After this showering phase the final-state 
consists of a set of partons with constituent masses. For processes involving the production 
and decay of unstable particles, including decay chains, rather than attempting to calculate 
high multiplicity matrix elements, the simulation is simplified by appealing to the narrow width 
approximation, i.e. treating the production and decay processes according to separate matrix 
elements, assuming no interference between the two. Unstable coloured particles are therefore 
produced in hard processes and the decays of other unstable particles, and showered like any other 
final-state coloured particle. In this case the showering process does not assign a constituent mass 
to the final state of the shower, but rather preserves whatever mass was assigned at the production 
stage. 
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For very high mass coloured particles, e.g. the top quark, the phase space available for the de- 
cay can be so large that the decay occurs before any hadronization can take place. Consequently, 
as well as undergoing time- like showers (q 2 > m 2 ) in their production phase, these partons will 
also undergo a further space-like showering (q 2 < m 2 ) of QCD radiation prior to their decay. In 
addition, due to colour conservation, the decay products themselves will also give rise to time-like 
showers. 

Since, in the narrow width approximation, the matrix element factorizes into one for the pro- 
duction process and another for the decay process, we may regard these as two independent hard 
processes, and this is the sense in which we simulate the associated parton showers. Given this 
picture it is immediately clear that the time-like parton showers, from coloured decay products, 
have an identical evolution to those used to simulate final-state radiation in the production pro- 
cess. Only the initial conditions for the shower evolution are different, although their selection 
is, nevertheless, still based on examining the colour flow in the underlying hard decay process 
(see Sect. 6.3.4). 

Conversely, the initial-state space-like shower created by a decaying particle is quite different 
to that of an initial-state particle from the production process (Sect. 6.5). In particular, it 
involves no PDFs, since the heavy parton originates from a hard scattering as opposed to a 
hadron. Furthermore, in the hard process it was necessary to evolve the initial-state partons 
backwards from the hard scattering to the incident hadrons, to efficiently sample any resonant 
structure in the underlying matrix elements. On the contrary, in decay processes, degrading 
the invariant mass of the decaying particle, via the emission of radiation, does not affect the 
efficiency with which any resonant structures in the decay matrix element are sampled. Hence, 
it is natural for the evolution of space-like decay showers to start with the unstable particle from 
the production process, and evolve it forward, towards its decay. 



6.6.1 Evolution 

As in our discussion of the other showering algorithms, the description here uses the Sudakov 
decomposition of the momenta given in Eq. (6.1). In space-like decay showers, the decaying 
particle ij undergoes branchings ij — * i + j, where j is a final-state time-like parton and % is the 
same decaying particle with an increased space-like virtuality: qf < q~. < m~.. In this process 
the original particle acquires a space-like virtuality, 



£ = z<% + P \^f, (6.85) 

where z = cxi/a-, p\ = —p 2 ± > 0, and p± = q±i — zc L±_ij- Since, in the decay shower, rrii = m^, 
the space-like evolution variable in Eq. (6.8) simplifies to 

„ „ zm 2 -p\ 
q 2 = m 2 + 1 F± . (6.86) 
{l-z) 

Unlike the previous discussions of final- and initial- state showers, here, by evolving forward 
toward the decay process, the evolution variable is increasing. The requirement that the relative 
transverse momentum of the branching is real, > 0, imposes an upper limit, z +) on z where 

9 

z + = 1 + 2 ( ~ 2 _ J m?) (l ~ y/l + 4 (<f - m 2 ) /mfj ■ (6.87) 
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For the space-like decay shower we have the further constraint that the parton showering cannot 
degrade the invariant mass of the decaying object below the threshold required for the decay 
process, which imposes a lower limit on z. 

Since no PDF is involved in this forward parton-shower evolution algorithm, the Sudakov 
form factor has exactly the same form as that used for final-state radiation in Eqs. (6.20,6.21). 
Consequently the forward evolution can be performed using the veto algorithm in almost exactly 
the same way as was done for the final-state showers (Sect. 6.4.1). The main difference is in 
the implementation of the angular ordering bounds for subsequent branchings. For final-state 
radiation involving branchings ij — * % + j, where i has a light-cone momentum fraction z, we 
evolved i and j downward from <jhi = zq and q~hj = (1 — z) q respectively, where q was the scale 
of the ij branching. Since the decay shower is really a forward-evolving initial-state shower, we 
evolve i upward from qh% = q and j downward from q^j = (1 — z) q. This procedure is iterated 
until the scale q approaches the minimum compatible with the threshold for the underlying decay 
process. 

6.6.2 Kinematic reconstruction 

In the approach of [20], for the simulation of QCD radiation in particle decays, the recoil due to 
the radiation emitted from the decaying particle is absorbed by its final-state colour partner. The 
reconstruction described in [23], valid for the decay of a coloured particle to a colour connected 
final-state particle and a colour-singlet system, was designed to preserve the mass of the colour- 
singlet system. In the case of top decay this amounts to preserving the mass of the W boson, 
and the momentum of the decaying particle. More complicated colour structures, involving more 
coloured particles in the final-state, e.g. gluino decays, require a generalization of this momentum 
reconstruction procedure. 

Consider the decay of a coloured particle with momentum p, to n + 1 particles. We denote the 
momentum of the colour partner of the decaying particle p, and the momenta of the remaining 
primary decay products are denoted Pi=i, n - Prior to simulating the effects of QCD radiation, 

p = p + J2Pi- (6-88) 

i=l 

After simulating parton-shower radiation in the decay, the original momenta of the decay products 
must be shifted and rescaled to accommodate the additional initial- state radiation. We require 
the sum of the new momenta of the colour partner, q, the other primary decay products, qi, and 
the radiation emitted prior to the decay, qisn, to equal that of the decaying particle: 

n 
i=l 

To achieve this momentum balance we rescale the three-momenta of all Pi by a common factor 
k\, and the three- momentum of the colour partner p by a separate factor hi. The component 
of the momentum of the emitted radiation transverse to the colour partner is absorbed by the 
colour partner. In the rest frame of the decaying particle these rescalings and shiftings look as 
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follows: 



p = (0;m); (6.90a) 
ft = (hPuJkflPtf+pi) ; (6.90b) 



9 = (MjP - q±/5i?; \/ k 2 |p| 2 + |q±/si?| 2 + P 2 J , (6.90c) 

where m is the mass of the decaying particle and q±isn is the component of the three-momentum 
of the initial-state radiation perpendicular to p. 

The rescaling factors k\ t 2 allow for the remaining conservation of energy and of momentum 
in the longitudinal direction. Three-momentum conservation in the longitudinal, p, direction 
requires that 

n 

hp + h Pi + <1\\isr = 0. (6.91) 

i=l 

The momentum of the initial-state radiation perpendicular to the direction of the colour partner, 
<Ilisr, can be projected out, leaving the parallel component <1\\isr, by taking the dot product 
with the spatial component of the n basis vector (aligned with p). Doing so gives 

h = k 2 + (6.92) 
p ■ n 

Finally, from the conservation of energy we have 



\Ail p *| 2 +p t + V k * 'P |2 + \^sr\ 2 +P 2 + E lSR = m, (6.93) 



i=l 



where -Eisr is the energy of the initial-state radiation. This system of equations Eqs. (6.91, 6.92, 
6.93) for the rescaling factors can be solved analytically for two-body decays, or numerically, 
using the Newton-Raphson method, for higher multiplicities. 



6.7 The running coupling constant as 

The running coupling constant enters every dynamical aspect of the parton shower, so a thorough 
treatment of it is mandatory for all parton shower simulations. 



6.7.1 The argument of q s 

As was noted in Sect. 6.2, our definition of the momentum fraction z is consistent with that used 
in the derivation of the quasi-collinear splitting functions, hence n does not just define a basis 
vector in the Sudakov decomposition but it also specifies the choice of light-cone (axial) gauge. 

Axial gauges have many special properties, most notable of these is that they are ghost-free. 
Another, related, interesting feature of the light-cone gauge is that, similar to QED, where the 
Ward identities guarantee the equality of the electron field and vertex renormalization constants, 
in light-cone gauge QCD, the Ward identities reveal that the 3-gluon vertex renormalization 
constant Z^3, is equal to that of the transverse components gluon field Z]( 2 [71]. This simplifies 
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the usual relation between the bare coupling g s and renormalized coupling constant gs from 
g$ = Z^3Z A 3 ^ 2 gs, to gf^ = Z~± gs, ^.e. in the light-cone gauge, the running of the QCD 
coupling constant is due to the gluon self-energy corrections alone. It is therefore no surprise 
that explicit, dimensionally regulated, one- loop calculations of the gluon self-energy in this gauge 
possess an ultraviolet divergence proportional to the usual QCD beta function [71,72]. 

In calculating higher order corrections to the splitting functions one must consider self-energy 
corrections to the emitted gluons and their associated counter-terms. The self-energy corrections 
are equal to zero because the gluons are on-shell and so the associated loop integrals have no 
scale, which means they vanish in dimensional regularization. This vanishing is essentially a 
complete cancellation of the ultraviolet and infrared parts of the integrals. Therefore including 
the counter-terms cancels explicitly the ultraviolet divergent parts of the loop integrals leaving 
behind infrared divergent parts, which must have the same pole structure as the ultraviolet parts 
i.e. they must also be proportional to the beta function. The residual virtual infrared divergence 
is canceled by the associated real emission corrections, in this case the two graphs where the 
emitted gluon splits either to two on-shell gluons or to a quark- ant iquark pair. 

As usual, this cancellation of infrared poles generates an associated logarithm, with the same 
coefficient as the pole (the beta function), of the phase space boundary divided by \i (the renor- 
malization scale) [54,73]. The phase space boundary is equal to the maximum possible virtuality 
of the daughter gluon, the branchings of which comprise the real emission corrections. For a 
time-like splitting, ij — > i + j where ij is a quark, % is a daughter quark and j is the daughter 
gluon, to which we consider real and virtual corrections, a quick calculation in the Sudakov basis 
Eq. (6.1) shows 

q 2 j <{l-z)q} j . (6.94) 

The net effect of these real and virtual corrections is therefore to simply correct the leading order 
q —> qg splitting function by a multiplicative factor 

1 - (3 a s {fi 2 ) In ((1 - z) q}J /) + O (a s ) , (6.95) 

where the omitted O (as) terms are non-logarithmic, non-kinematic, constant terms, (3$ is the 
QCD beta function, and /i 2 is the renormalization scale. 

Two important points follow directly from this analysis. Firstly, for soft configurations, z — > 1, 
the effect of these loop contributions can produce large, numerically significant, logarithms. 
Secondly, plainly, by choosing the have 

q® = a®P® + (3®P® + q±®- (6.96) 

renormalization scale to be (1 — z) o~, instead of the more obvious g~, the corrections vanish, or 
rather, more correctly, they are absorbed in the coupling constant. 

For g — > gg splittings the same arguments hold but in this case it is apparent that as well 
as large logarithms of 1 — z, large logarithms of z are also possible from soft emission in the 
z ^ region. We may simultaneously include both types of correction by using z (1 — z) q~. as 
the argument of the running coupling, which we implement in practice as 

a s {z 2 {l-zfq 2 ). (6.97) 

From the point of view of the leading-log approximation, the choice of scale is technically a 
higher order consideration, nevertheless, these effects turn out to be highly phenomenologically 
significant, particularly their effect on multiplicity distributions and cluster mass spectra [73,74]. 
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6.7.2 The Monte Carlo scheme for a s 



We reiterate that by choosing the scale of the running coupling as advocated in Sect. 6.7.1 
(Eqs. (6.94,6.97)) we have 

hma s [(l-z)qf j jPj 1 l qg (z) = a s ^ (1 - a s p ln(l - z)) + 0{a 3 s ), (6.98) 

where we have momentarily abbreviated as ylfjj by and used a superscript [1] to denote 
that P [ ql qg is the one-loop [i.e. leading order) q — > qg splitting function. This is almost, but not 

[21 

exactly equal to the soft z — > 1 singular limit of the two-loop q — > qg splitting function PqX qg 

with a q evaluated at o~, 

y 

lim a s Pfi qg (z) = as ^(l- a s (3 In (1 - z) + ^-K g ) + <D(a s ), (6.99) 
where 19 

K * = c *(%-T)- T * n 'lF- (6 ' 100) 

On integrating over the phase space of the two-loop splitting function the K g term gives rise 
to terms ~ a'g\n 2 q~, i.e. it gives next-to- leading log soft-collinear contributions to the Sudakov 

exponent ~ aghi n q~. (as opposed to leading- log contributions ~ agln n+1 q~.). In a similar way 
to that in which the higher order ftoas In (1 — z) term was included, we may exploit the fact that 
the z — > 1 dependence of the K g term in PfX qg (z) is equal to that of PfX qg {z), to incorporate it 
in the running coupling as well. 

This is done by swapping the usual Aj^g QCD scale, from which the coupling runs, for 
A MC [64], 

Amc = Ams exp (K g /4n/3 ) , (6.101) 

where MC denotes the so-called Monte Carlo scheme. Expanding asP^Xqg (z) again, as in 
Eq. (6.98), but with as evaluated at (1 — z) q~. in the MC scheme, reproduces exactly the two- 
loop result in Eq. (6.99). With this prescription the Sudakov form factor generally includes all 
leading and next-to-leading log contributions, except for those due to soft wide angle gluon emis- 
sions, however, for the case that the underlying hard process comprises of a single colour dipole, 
these are also included (see Sect. 6.2 and [32,65]). 



6.7.3 Options for the treatment of 0:5 in parton showers 

Although we have made strong physical arguments restricting the argument of the coupling 
constant and suggesting a more physical renormalization scheme, there is still some degree of 
freedom in how precisely as is calculated. In what follows below we enumerate the options 
associated with these in the program. 

19 In fact the constants K g are given by the finite remainder of the real emission phase space corrections due to 
the daughter gluon splitting discussed in the last Sect. 6.7.1 (see e.g. Eqs. (5.28,C.12,C.13) of [75]). 
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InputOption This option selects the way in which initial conditions for running the coupling 
constant are determined. The default setting [InputOption=AlphaMZ] uses the experimen- 
tally determined value of as at the Z° resonance to calculate a value of Aqcd from which 
to run the coupling constant. This experimental input can be reset from the default value 20 
of 0.127 using the AlphaMZ interface. Alternatively one may select an option 
[InputOption=LambdaQCD] , which uses the input or default value of A^jg to calculate the 
coupling. The default value used for Aj^g- is 0.208 GeV, which may be reset using the interface 
LambdaQCD. 

LambdaOption This option determines whether the value of Aqcd ; calculated from «5 (m^o) 
or input according to InputOption, is given in the MC (Monte Carlo) scheme of Ref. [64], as de- 
scribed in Sect. 6.7.2 [LambdaOption=Same] , the default, or the MS scheme 
[LambdaOption=Convert] . 

NumberOfLoops This parameter specifies the loop order of the beta function used to calculate 
the running of as- The default setting uses the three-loop beta function. 

ThresholdOption This option selects whether to use the current 
[ThresholdOption=Current] or constituent [ThresholdOption=Constituent] quark 
masses in determining the flavour thresholds in the evolution of the coupling constant. The 
default setting uses the (MS) current quark masses. 

Qmin The Qmin parameter represents the scale beneath which non-perturbative effects are 
considered to render the usual renormalization group running with a beta function determined 
at some finite loop order, invalid. Below this scale, which is currently tuned to 0.935 GeV, a 
number of parameterizations of the scaling of the coupling with energy may be selected according 
to the NPAlphaS option described below. 

NPAlphaS The NPAlphaS option selects a parameterization of the scaling of the running 
coupling with energy in what we regard as the non-perturbative region, where the scale at 
which it is to be evaluated falls below the value set by Qmin. By setting [NPAlphaS=Zero] 
the coupling is simply taken to be zero for scales Q < Qmin. For [NPAlphaS=Const] the 
coupling freezes out at Qmin, i.e. it assumes the constant value as = as (Qmin) for all scales 
below Qmin. This is the default parameterization. It is the same prescription used in early 
works on resummation by Curci and Greco [76,77]. The options [NPAlphaS=Linear] and 
[NPAlphaS=Quadratic] calculate the running coupling below Qmin according to asQ/Qmin 
and as (Q/Qmin) 2 respectively. Setting [NPAlphaS=Exxl] assumes a quadratically decreasing 
running of the coupling in the non-perturbative region from the value AlphaMaxNP down to 
as- Finally, [NPAlphaS=Exx2] sets as equal to AlphaMaxNP for all input scales Q < Qmin, 
which amounts to a minor variation of the default freeze-out option. 

20 The default value is tuned to e + e~ annihilation data as described in Appendix D and is typical of the values 
one gets when fitting leading order QCD predictions to data. 
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6.8 Matrix element corrections 



As stated in Sect. 6.2, the effects of unresolvable gluon emissions have been included to all orders 
through the Sudakov form factor. The master formula and shower algorithms generate further 
resolvable emissions by approximating the full next-to-leading order real emission matrix element 
by a product of quasi-collinear splitting functions multiplying the tree level amplitude. Ideally, 
we wish to include higher-order effects as accurately as possible and do this for certain processes 
using matrix element corrections. We aim to correct two deficiencies of the shower algorithm: 
(i) it may not cover the whole phase space, leaving a region of high p± (i.e. non-soft non-collinear) 
emission unpopulated; and (ii) even in the region it does populate, as one extrapolates away from 
the soft and collinear limits it may not do a perfect job. We call these the hard and soft matrix 
element corrections respectively [78]. 



6.8.1 Soft matrix element corrections 

In the parton shower approximation the probability density that the zth resolvable parton is 
emitted into [q 2 , q 2 + dq 2 ] , [z, z + dz] is 

dV & ^ =^T dZ P ^ ij ( "' f) 9 (P1 " 0) • (6 ' 102) 

This approximation works well for the case that the emission lies within the domain of the 
quasi-collinear limit. On the other hand the exact matrix element calculation gives us that the 
probability of a resolved emission as 

f dV m c - = [ dq 2 dz 6 (pi > 0) , (6.103) 

Jn J o-Qdzdq 2 

where da is the differential cross section for the underlying process with a further parton emis- 
sion, and 1Z denotes the region of phase space corresponding to resolved emissions. The KLN 
and Bloch-Nordsieck theorems imply that all large logarithmic corrections to the cross section 
must vanish once the full available phase space is integrated over. It follows that the O (as) 
contribution to the total cross section from an unresolved emission may be written — j n dP m,c ', 
at the level of large (leading and next-to-leading) logarithms. Proceeding in the same way as our 
earlier derivations Eq. (6.1), we then have that the probability density that the zth resolvable 
gluon is emitted into [q 2 , q 2 + dq 2 ], [z, z + dz] is given by the integrand of 

^-i 1 d 2 a ( 1 d 2 a 

d $ dz — T~^2 exp - / d^de — — — I . (6.104) 

2 . cr dzdqf \ J;2 a dzdq J 



We may generate the distribution in Eq. (6.104) by simply augmenting the veto algorithm that 
is used to produce Eq. (6.21) with a single additional rejection weight, simply vetoing emissions 
if a random number IZs is such that 



dP 1 

n s > — 
s - dV 



(6.105) 



For this to work we require that the parton shower emission probability dP always overestimates 
that of the exact matrix element dP m e ', if necessary this can be achieved by simply enhancing 
the emission probability of the parton shower with a constant factor. 
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This correction is consistently applied to every emission that has the highest so far in the 
shower. This ensures not only that the leading order expansion of the shower distribution agrees 
with the leading order matrix element, but also that the hardest (i.e. furthest from the soft and 
collinear limits) emission reproduces it. One might be concerned that it is really only proper 
to apply this correction to the final, largest p^emission, however, in the context of a coherent 
parton branching formalism (angular ordering) the earlier wide-angle emission is considered too 
soft to resolve the subsequent, smaller angle but larger splitting, and is therefore effectively 
distributed as if the latter emission did not occur. In this way, not only the hardest emission 
is improved by the correction, but all reasonably hard wide-angle emissions. Thus the correct 
procedure is to correct all those emissions that are the hardest so far, from the distribution in 
Eq. (6.21) to that in Eq. (6.104) by applying the veto in Eq. (6.105) [78]. 

Given that each soft matrix element correction amounts to exponentiating the next-to-leading 
order real emission matrix element divided by the leading order matrix element, provided one 
selects the option to evaluate the running coupling in the Monte Carlo scheme [64] , the Sudakov 
form factor is in this case formally of next-to-leading log accuracy for corrections to processes 
comprised of a single colour flow 21 . For processes involving more than one underlying colour the 
next-to-leading log accuracy of the Sudakov form factor is only correct up to terms O (l/N^) 
[32,65]. 



6.8.2 Hard matrix element corrections 

In addition to correcting the distribution of radiation inside the regions of phase space that are 
populated by the parton shower, we also wish to correct the distribution of radiation outside, in 
the high p^, unpopulated, dead region. We wish to distribute the radiation in the dead regions 
according to the exact tree-level real emission matrix element i.e. according to 

1 /*«£ i,max j .max ip^i) H^fT 

— / dxi / d Xj — - — , (6.106) 

(J() 1^. fry.. ■ (rf-\ (IXjdXj 

u " • ij i,min ' J • L, j,min\' ij i} " J 

where da is the differential cross section obtained using the next-to-leading order, real emission 
matrix element, and (xi,Xj) are variables parameterizing the phase space associated with the 
emission of the extra parton. 

The algorithm for populating the dead region is basic in principle. Prior to any showering the 
program checks if a matrix element correction is available for the hard process. If one is available 
the algorithm then generates a point in the appropriate region of phase space, ideally with some 
importance sampling of the integrand. The differential cross section associated with this point, 
as given in Eq. (6.106), is evaluated and multiplied by a phase space volume factor V (x,) given 
by 

V {p^i) i.^i,max ^i,mm) (^j,max (-^i) ■Ej,min j (6.107) 

giving the event weight. The emission is retained if this weight is less than a uniformly distributed 
random number 1Z G [0, 1], and the momenta of the new parton configuration are reconstructed 
from the generated values of and Xj . 



21 For processes involving initial-state radiation, this also requires evaluating the parton densities at a scale of 
order pj_ [32]. 
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6.8.3 Using HerwigH — \- matrix element corrections 

The current version of Herwig++ contains matrix element corrections for four different hard 
processes: neutral and charged current Drell-Yan processes, gg — > h°, top quark decays and 
e + e~ — > gg processes. The associated C++ classes are DrellYanMECorrection, GGtoHMECorrection, 
TopDecayMECorrection and VectorBosonQQbarMECorrection. 

Naturally each of these pro cess- dependent matrix element corrections checks whether it cor- 
responds to the hard process (or, for top quark decays, the decay process). In other words, users 
need not worry that, if matrix element corrections are globally switched on in the code, the 
correction for e.g. the Drell-Yan processes is applied to the gg — ► H process they have selected 
to generate. 

All three corrections are loaded in the Repository in the default set-up. The switch 
MECorrMode determines the way in which all matrix elements are used. If [MECorrMode=0] 
is selected no matrix element corrections will be applied at all. The default setting 
[MECorrMode=l], applies both the hard and soft matrix element corrections for each one 
loaded in the Repository (if the associated processes are generated). Options [MECorrMode=2] 
and [MECorrMode=3] turn off the soft and hard matrix element corrections respectively. 

6.9 Showering in the POWHEG scheme 

In the POWHEG approach [31] the NLO differential cross section for a given TV-body process can 
be written as 

dcr = B (<& B ) d$e 
where B ($b) is defined as 

B($ B ) =B^ b ) + V(^ b ) + J |£($b,$*)- X}Ci($B,$*) J d<f R , (6.109) 

B(<&b) is the leading-order contribution, dependent on the A^-body phase space variables 
the Born variables. The regularized virtual term V ($b) is a finite contribution arising from 
the combination of unresolvable real emission and virtual loop contributions. The remaining 
terms in square brackets are due to A+l-body real emission processes which depend on both the 
Born variables and additional radiative variables, parametrizing the emission of the extra 
parton. The real emission term, R($b, &r), is given by a sum of parton flux factors multiplied 
by real emission matrix elements for each channel contributing to the NLO cross section. Finally, 
each term Cj (<&b, $_r) corresponds to a combination of real counterterms / counter- event weights, 
regulating the singularities in R The modified Sudakov form factor is defined as 

A r (Pt) = exp 

where hr ($b, $r) is equal to the transverse momentum of the extra parton. 

As well as circumventing the problem of negative event weights the POWHEG method defines 
how the highest pt emission may be modified to include the logarithmically enhanced effects 
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(6.108) 
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(6.110) 
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of soft wide-angle radiation. In Ref. [31] it was shown how the angular-ordered parton shower 
which produces the hardest emission, may be decomposed into a truncated shower simulating 
coherent, soft wide-angle emissions, followed by the highest px (hardest) emission, followed again 
by further vetoed parton showers, comprising of lower px, smaller angle emissions. Performing 
this decomposition established the form of the truncated and vetoed showers, thereby describing 
all of the ingredients necessary to shower the radiative events in the POWHEG approach. This 
procedure was proven in [32] to give agreement with the NLO cross section, for all inclusive 
observables, while retaining the logarithmic accuracy of the shower. 

In the POWHEG framework positive weight events distributed with NLO accuracy can be 
showered to resum further logarithmically enhanced corrections by: 

• generating an event according to Eq. 6.108; 

• directly hadronizing non-radiating events; 

• mapping the radiative variables parametrizing the emission into the evolution scale, mo- 
mentum fraction and azimuthal angle (q^, Zh, <ph)i from which the parton shower would 
reconstruct identical momenta; 

• using the original leading-order configuration from B ($#) evolve the leg emitting the extra 
radiation from the default initial scale, determined by the colour structure of the iV-body 
process, down to the hardest emission scale q^ such that the px is less than that of the 
hardest emission px h , the radiation is angular-ordered and branchings do not change the 
flavour of the emitting parton; 

• inserting a branching with parameters (q~h, z^, (f>h) into the shower when the evolution scale 
reaches qh, 

• generating px vetoed showers from all external legs. 

This procedure allows the generation of the truncated shower with only a few changes to the 
normal Herwig++ shower algorithm. 

In the Herwig++ implementation the generation of the Born variables according to Eq. 6.108 
is performed as described in Sect. 3.2. The rest of the POWHEG algorithm is then implemented by 
using the PowhegEvolver class which inherits from the Evolver class and implements the generation 
of the hardest emission and truncated shower. 

The hardest (highest px) emission is generated from the iV-body configuration according to the 
modified Sudakov form factor, Eq. 6.110, using a class inheriting from HardestEmissionGenerator 
which implements the hardest emission for a specific process. Currently only gauge bosons 
via the Drell-Yan process, including virtual gauge bosons in processes like the production of a 
gauge boson in association with the Higgs boson, and Higgs production via gluon fusion are 
implemented. 

In order to perform the truncated shower this emission is then interpreted as an emission from 
the parton shower in the following way 22 . This is essentially the inverse of the reconstruction 
of the parton shower: first the reshuffling boosts are inverted; then the shower variables are 

22 We only describe the case of initial-state emission together with an initial-state colour partner as the only 
processes currently implemented are of this type. 
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calculated by decomposing the momenta after the boost in terms of the Sudakov basis. The 
momenta of the potentially off-shell partons which would have participated in the hard process 
are calculated, for initial-state radiation these can be decomposed in terms of the Sudakov basis 
as 

Q® = + P'®P® + <?!©• (6.111) 

We then need to calculate the on-shell momenta in the iV-body hard process which would generate 
them, and the shower variables for the hard emission. The momentum fractions and the partons 
before the emission, x@, can be calculated from the requirement that the centre-of-mass energy 
and rapidity of the collision are preserved. The rescaling parameters in Eq. 6.73 are 

fee = ^, k e = ^, (6.112) 

which allows the inverse of the boosts applied in the momentum reshuffling to be calculated. 
The momenta the radiation parton would have in the shower before the reshuffling can then 
be calculated by performing the inverse of the rescaling boost. The momentum fraction of the 
emission is given by 

z=—, (6.113) 

where on is the Sudakov parameter for the space-like parton entering the hard process and 
the Sudakov parameter of the initial-state parent parton. In this simple case the transverse 
momentum is simply equal to that of the off-shell space-like parton initiating the leading-order 
hard process, or equivalently, its outgoing, time-like, sister parton. The scale of the branching is 
defined in terms of the pt and light-cone momentum fraction z, as 

f - %±4. (6.H4) 
(1 - z) 

This procedure is implemented in the QTildeReconstructor class. 

Once this inverse reconstruction has been performed the PowhegEvolver can shower the event 
including both truncated and vetoed showers using information on the hardest emission stored 
in a HardTree object which is produced by the HardestEmissionGenerator. 



6.10 Code structure 

The Herwig++ shower module consists of a large number of classes and is designed to be flexible, 
in the sense that any DGLAP-type shower evolution based on 1 — > 2 branchings where momentum 
conservation is enforced globally after the evolution has been performed can be implemented. The 
only concrete implementation so far is the improved angular-ordered shower based on [20] and 
described above. 

We will only describe the structure of the code, i.e. how the various classes work together to 
generate the parton shower evolution. Detailed documentation of all the classes can be found in 
the Doxygen documentation. In a future release, the structure will be slightly changed to allow 
for more general shower evolution, such as dipole-type showers. 

The main class implementing the Herwig++ shower is the ShowerHandler class, which inherits 
from the CascadeHandler class of ThePEG. It has responsibility for the overall administration of 
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the multiple interactions, as described in Sect. 8, the showering of primary and secondary hard 
scattering processes, the decay of any unstable fundamental particles 23 and the generation of 
any radiation produced in their decays. The ShowerHandler uses a number of helper classes to 
implement various parts of the algorithm together with some data storage classes, which hold 
information needed to generate the parton shower. 
The ShowerHandler proceeds as follows: 

• The Event object supplied to the ShowerHandler is first analysed and the particles to be 
showered extracted. These particles are converted from Particle objects, which store particle 
information in ThePEG, to ShowerParticle objects, which inherit from Particle and include 
the storage of the additional information, such as the evolution scales and colour part- 
ners, needed to generate the parton shower. Each particle in a hard process, be that the 
primary scattering process or the subsequent decay of a fundamental particle, is assigned 
to a ShowerProgenitor object containing references to the particle together with additional 
information required for particles that initiate a parton shower. For each hard process a 
ShowerTree object is created containing the ShowerProgenitor objects for all the particles in 
the hard process and the information required to shower that process. 

• The ShowerHandler uses the helper Evolver to generate the radiation from each hard scat- 
tering or decay process. Once the parton showers have been generated for all the hard 
processes the ShowerHandler inserts them into the Event object. 

• The MPIHandler then generates any secondary hard scatterings required, which are subse- 
quently showered by the Evolver, as described in Sect. 8. 

• Finally, after all the scatterings have been showered, the hadronic remnant is decayed to 
conserve momentum and flavour using the HwRemDecayer class. 

The main helper class of the ShowerHandler is the Evolver, which is responsible for generating the 
parton shower from an individual hard process, stored as a ShowerTree object. The Evolver first 
finds the colour partners and initial scale for the parton showers from each particle, as described 
in Sect. 6.3. At this stage, if there is a suitable class inheriting from MECorrectionBase, which 
implements the matrix element correction for the process as described in Sect. 6.8, the hard 
matrix element correction is applied. The Evolver is also currently responsible for generating the 
intrinsic pj_ of incoming partons in hadronic collisions at this stage. 

Given the initial scale, the evolution of the particles proceeds as described in Sects. 6.4-6.6, 
using the SplittingGenerator class to generate the types and scales of the branchings. In turn the 
SplittingGenerator uses the SudakovFormFactor to generate the possible evolution scales for each 
allowed type of branching and then selects the branching with the highest scale, as described in 
Sect. 6.4. The new ShowerParticles produced in the branching are then evolved until no further 
branching is possible. When all the particles have been evolved the KinematicsReconstructor 
reconstructs the momentum of all the particles in the shower (Sects. 6.4-6.6). 

The ShowerHandler and Evolver classes are mainly administrative, the actual physics is im- 
plemented in the various helper classes. For this reason these helper classes, which are specific 

23 Currently most fundamental particle decays are performed before the parton shower is generated, although 
in future we plan to generate them as part of the parton-shower algorithm. 
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to the details of the parton shower algorithm, are contained in the ShowerModel class. It is in- 
tended that different DGLAP based parton shower algorithms, for example the original angular- 
ordered parton shower algorithm used in FORTRAN HERWIG, can be implemented by inheriting 
from the ShowerModel and specifying the helper classes to be used in that model, which in- 
herit from the KinematicsReconstructor, PartnerFinder, SudakovFormFactor and MECorrectionBase 
classes. For example, the QTildeModel, which implements the improved angular-ordered shower 
described above, uses the QTildeReconstructor, QTildeFinder, QTildeSudakov and 
QTildeMECorrection classes. 

In turn many of the helper classes used by the main classes implementing the shower have 
their own helper classes for various parts of the simulation. 

The SplittingGenerator class holds lists of available branchings, providing interface switches to 
either enable or disable radiation, in the initial or final state, for different interactions. They are 
used to generate the shower variables associated with each branching using SudakovFormFactor 
objects. The SplittingGenerator and SudakovFormFactor classes use the following helper classes: 

SplittingFunction This is the base class for defining splitting functions used in the shower 
evolution. This includes the calculation of the splitting function together with the overestimate, 
integral and inverse integral of it required to implement the veto algorithm as described in 
Sects. 6.4 and 6.5. The splitting functions implemented in Herwig++ are listed in Sect. 6.2. 

ShowerAlpha This is the base class implementing the running couplings used in the shower 
evolution. 

The Evolver uses the ShowerVeto class to provide a general interface to veto emission attempts 
by the shower. The veto may be applied to either a single emission (resetting the evolution scale 
for the particle to the attempted branching scale), an attempt to shower a given event, or the 
overall event generation. 

The additional features needed in the POWHEG approach are implemented in the 
PowhegEvolver class which inherits from the Evolver class and uses classes inheriting from 
HardestEmissionGenerator to generate the hardest emission in the POWHEG scheme. 

Finally three special exception classes are used inside the shower module, mainly to communi- 
cate exceptional events or configurations, rather than signaling a serious error during event gen- 
eration. The exceptions are handled completely within the shower module. In particular we use 
VetoShower to communicate vetoing of a complete shower attempt. 
KinematicsReconstructionVeto is used to signal an exceptional configuration that cannot be han- 
dled by the KinematicsReconstructor, resulting in restarting the shower from the original event 
(similar to a VetoShower exception). ShowerTriesVeto signals that complete showering of a given 
event failed a predefined number of times. This is handled together with the generation of 
multiple interactions. 
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7 Hadronization 



After the parton shower, the quarks and gluons must be formed into the observed hadrons. The 
colour preconfmement property [74] of the angular-ordered parton shower is used as the basis of 
the cluster model [2], which is used in Herwig++ to model the hadronization. This model has the 
properties that it is local in the colour of the partons and independent of both the hard process 
and centre-of-mass energy of the collision [2,3]. 

7.1 Gluon splitting and cluster formation 

The first step of the cluster hadronization model is to non-perturbatively split the gluons left 
at the end of the parton shower into quark- ant iquark pairs. Since, at the end of the Herwig++ 
shower the gluons are given their constituent mass it is essential that this mass is heavier than 
twice the constituent mass of the lightest quark 24 . The gluon is allowed to decay into any of the 
accessible quark flavours with probability given by the available phase space for the decay 25 . 

The gluon decays isotropically and following this isotropic decay the event only contains colour 
connected (di) quarks and anti-(di) quarks. The colour singlets formed by these colour connected 
parton pairs are formed into clusters with the momentum given by the sum of the momenta of 
the constituent partons. The principle of colour-preconfinement states that the mass distribution 
of these clusters is independent of the hard scattering process and its centre-of-mass energy. As 
can be seen in Fig. 4a, the shower algorithm in Herwig++ obeys preconfmement fairly well by 
100 GeV and is clearly invariant beyond that. 

7.2 Cluster fission 

The cluster model is based on the observation that because the cluster mass spectrum is both 
universal and peaked at low masses, as shown in Fig. 4a, the clusters can be regarded as highly 
excited hadron resonances and decayed, according to phase space, into the observed hadrons. 
There is however a small fraction of clusters that are too heavy for this to be a reasonable 
approach. These heavy clusters are therefore first split into lighter clusters before they decay. 
A cluster is split into two clusters if the mass, M, is such that 

M cip OW > Cl max clpow + (mi + m 2 ) clpow , (7.1) 

where Cl max and Cl pow are parameters of the model, and mi j2 are the masses of the constituent 
partons of the cluster. In practice, in the most recent version of the model, in order to improve 
the description of the production of bottom and charm hadrons, we include separate values of 
both Cl max (CIMaxLight, CIMaxCharm and CIMaxBottom) and Cl pow (CIPowLight, 
CIPowCharm, CIPowBottom) for clusters containing light, charm and bottom quarks re- 
spectively. The default values of these and other important hadronization parameters are given 
in Table 10 at the end of this Section. 

24 We normally take the constituent masses of the up and down quarks to be equal although they can in principle 
be different. 

25 The option of gluon decay into diquarks, which was available in FORTRAN HERWIG, is no longer supported. 
Diquarks are therefore present only as remnants of incoming baryons, or from baryon number violating pro- 
cesses (see Sect. 7.4.2). 
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Figure 4: The mass spectrum of a) the primary clusters and b) the clusters after cluster fission. 
The solid, dashed and dot-dashed lines show the clusters produced in hadronization of e + e~ — > dd 
events at a centre-of-mass energy of 100 GeV, 1 TeV and 10 TeV respectively. Only clusters 
containing light quarks are shown. 



For clusters that need to be split, a qq pair is selected to be popped from the vacuum. Only 
up, down and strange quarks are chosen with probabilities given by the parameters Pwtj 26 , 
where % is the flavour of the quark. Once a pair is selected the cluster is decayed into two new 
clusters with one of the original partons in each cluster. Unless one of the partons is a remnant 
of the incoming beam particle the mass distribution of the new clusters is given by 

Mi = m 1 + (M - mi - m q )n\ /P , (7.2a) 
M 2 = m 2 + (M — m 2 — m g )Tl l 2 /P , (7.2b) 

where m q is the mass of the parton popped from the vacuum and M 12 are the masses of the 
clusters formed by the splitting. The distribution of the masses of the clusters is controlled by the 
parameter P, which is PSplitLight, PSplitCharm or PSplitBottom for clusters containing 
light, charm or bottom quarks. 

In addition to the selection of the mass according to Eq. (7.2) the masses of the daughter 
clusters are required to be less than that of the parent cluster and greater than the sum of the 
masses of their constituent partons. The spectrum of the cluster masses after the cluster splitting 
is shown in Fig. 4b. 

For clusters that contain a remnant of the beam particle in hadronic collisions a soft distribu- 
tion is used for the masses of the clusters produced in the splitting. The RemnantOption switch 
controls whether the soft distribution is used for both daughter clusters [RemnantOption=0] 

26 Wc use Pwti to denote the probability of selecting a given quark or diquark. This is given by the parameters 
PwtDquark, PwtUquark, PwtSquark. PwtCquark and PwtBquark for the quarks and the product of the 
diquark probability PwtDiquark, the probabilities of the quarks forming the diquark, and a symmetry factor 
for diquarks. 
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or only the daughter cluster containing the remnant [RemnantOption=l], the default. The 
mass of the soft clusters is given by 

Mi = rrii + m q + x, (7.3) 
where x is distributed between and M — rri\ — m 2 — 2m q according to 

dP 

— = exp (-bx) , (7.4) 
dar 

where b = 2/SoftClusterFactor. 
7.3 Cluster decays 

The final step of the cluster hadronization model is the decay of the cluster into a pair of 
hadrons. For a cluster of a given flavour (91,92) a quark-antiquark or diquark-antidiquark pair 
(9, 9) is extracted from the vacuum and a pair of hadrons with flavours (91, q) and (9, 92) formed. 
The hadrons are selected from all the possible hadrons with the appropriate flavour based on the 
available phase space, spin and flavour of the hadrons. While the general approach is the same 
in all cluster models there are some variations. In Herwig++ the original model of Ref. [2] used 
in FORTRAN HERWIG [5,6], the approach of Ref. [79], which was designed to solve the problem 
of isospin violation in the original model if incomplete SU(2) multiplets of hadrons are included, 
and a new variant that addresses the issue of the low rate of baryon production in the approach 
of Ref. [79] , are implemented. 

In all these approaches the weight for the production of the hadrons and fe( 5) g 2 ) is 

W( fl («l,3)> b (Q,Q2)) = PqW a S a W b S b p* a>b , (7.5) 

where P q is the weight for the production of the given quark-antiquark or diquark-antidiquark 
pair, w a ,b are the weights for the production of individual hadrons and s a , b are the suppression 
factors for the hadrons, which allow the production rates of individual meson multiplets, and 
singlet and decuplet baryons to be adjusted. The momentum of the hadrons in the rest frame of 
the decaying cluster, 

Pl, b = 2^ [( M2 - K + m b ) 2 ) (M 2 - (m a - m b f)] * , (7.6) 

measures the phase space available for two-body decay. If the masses of the decay products 
are greater than the mass of the cluster then the momentum is set to zero. The weight for the 
individual hadron is 

w h = w mix (2J h + 1), (7.7) 

where w m i x is the weight for the mixing of the neutral light mesons 27 and Jh is the spin of the 
hadron. 

The different approaches vary in how they implement the selection of the cluster decay prod- 
ucts based on this probability. 



27 



w m ix = 1 for all other particles. 
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In the approach of Ref. [2] the probability is generated in a number of pieces. First the flavour 
of the quark- ant iquark, or diquark-antidiquark, pair popped from the vacuum is selected with 
probability 

Both the hadrons produced in the cluster decay are then selected from the available hadrons 
of the appropriate flavours using the weight 

(7.9) 



^max(g,g') 

where w max r qq n is the maximum value of the weight for a given flavour combination. 
A weight is calculated for this pair of hadrons 

W = !^k y ( 7.io) 

Pm&x 

where p* a b is the momentum of the hadrons in the cluster rest frame and p max is the maximum 
momenta of the decay products for hadrons with the relevant flavour 28 . The hadrons produced 
are then selected according to this weight. 

This procedure approximately gives a probability 

n a (*i,g)> 6 ( 9 ,fc)l9i>92) = T7 ~* I 7 - 11 ) 

(<?ii<?) iv (g,92) U7 max((ji ,q) U7 max(<j,<j2) Fmax 

of choosing hadrons cb{ qi ,q) an d ^(g,g 2 )- The number of hadrons with flavour (g 1; q 2 ) is N( qi>q - 2 y 

Kupco [79] pointed out one problem with this approach: as new hadrons with a given flavour 
are added, the production of the existing hadrons with the same flavour is suppressed. In order 
to rectify this problem he proposed a new approach for choosing the decay products of the 
cluster. Instead of splitting the probability into separate parts, as in Ref. [2], a single weight was 
calculated for each combination of decay products 

^( a (ffi,5)> 6 («,«)l?i>fc) = p q w a w b s a s b p* ab , (7.12) 
which gives the probability of selecting the combination 

P(a {quq - h b {M2) \ qi ,q 2 ) = — — — — -. (7.13) 

The addition of new hadrons now increases the probability of choosing a particular flavour, 
however because these new hadrons are usually heavy they will not contribute for the majority 
of light clusters. 

The main problem with this approach is that because many more mesons are included in the 
simulation than baryons not enough baryons are produced. In order to address this problem 
in Herwig++, if a cluster mass is sufficiently large that it can decay into the lightest baryon- 
antibaryon pair the parameter Pwt qq is used to decide whether to select a mesonic or baryonic 



2X 



That is, the momentum with the lightest possible choices for a and b. 
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decay of the cluster. The probabilities of selecting a mesonic decay or baryonic decay are 1+P 1 wt 

and jjp^f ■ This modification not only increases the number of baryons produced but gives 
direct control over the rate of baryon production. 

Once the decay products of the cluster are selected, the cluster is decayed. In general the 
cluster decay products are isotropically distributed in the cluster rest frame. However, hadrons 
that contain a parton produced in the perturbative stage of the event retain the direction of 
the parton in the cluster rest frame, apart from a possible Gaussian smearing of the direction. 
This is controlled by the CIDir parameter, which by default [ClDir=true] retains the parton 
direction, and the CISmr parameter, which controls the Gaussian smearing through an angle 
#smcar where 

cos# smoar = 1 + CISmr log 7?.. (7-14) 

The azimuthal angle relative to the parton direction is distributed uniformly. To provide greater 
control the parameters CIDir (CIDirLight, CIDirCharm and CIDirBottom) and CISmr 
(CISmrLight, CISmrCharm and CISmrBottom) can be set independently for clusters con- 
taining light, charm and bottom quarks. 

In practice there is always a small fraction of clusters that are too light to decay into two 
hadrons. These clusters are therefore decayed to the lightest hadron, with the appropriate 
flavours, together with a small reshuffling of energy and momentum with the neighbouring clus- 
ters to allow the hadron to be given the correct physical mass. The cluster with the smallest 
space-time distance that can absorb the recoil is used. In addition, for clusters containing a 
bottom or charm quark in order to improve the behaviour at the threshold the option exists of 
allowing clusters above the threshold mass, M threshold , for the production of two hadrons to decay 
into a single hadron such that a single hadron can be formed for masses 

M < M limit = (1 + SingleHadronLimit)M thrcshold . (7.15) 

The probability of such a single-meson cluster decay is assumed to decrease linearly for M thresho i d < 
M < Miimit. The parameters SingleHadronLimitCharm and SingleHadronLimitBottom 
control the limit on the production of single clusters for charm and bottom clusters respectively. 
Increasing the limit has the effect of hardening the momentum spectrum of the heavy mesons. 

7.3.1 Mixing weights 

For neutral mesons that only contain the light (up, down and strange) quarks there is mixing. If 
we consider the wavefunctions of the neutral mesons, which we write for the 1 Sq meson multiplet 
but the treatment applies to an arbitrary SU(3) flavour multiplet, then 

7T° = —=(dd — uu), (7.16a) 
v2 

r) = ipg cos 9 — tpi sin 9, (7.16b) 

r]' = ip 8 sin 9 + ip\ cos 9, (7.16c) 
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where 6 is the nonet mixing angle and the wavefunctions for the octet and singlet components 
are 



-08 = —= [uu + dd — 2ssj , (7.17a) 
V6 

ipi = —j= [uu + dd + ss) . (7.17b) 
v3 

The probabilities of finding a given quark-antiquark inside a particular neutral meson can be 
calculated, which gives the mixing weights for the neutral light mesons 

wf s = 0, (7.18a) 
w v sS = sin 2 (fl + 0), (7.18b) 
w ]' s = cos 2 (# + 0), (7.18c) 

where = tan -1 \/2 is the ideal mixing angle. 

In the approach of Ref. [2] the factor of | in the weights for the uu and dd components was 
omitted as this is approximately given by the ratio of the number of charged mesons containing 
up and down quarks to neutral ones, which is exactly two for ideal mixing where the ss mesons 
do not mix with those containing up and down quarks. 

In practice the mixing angles can be adjusted for each meson multiplet that is included in 
the simulation although with the exception of the lightest pseudoscalar, vector, tensor and spin-3 
multiplets the assumption of ideal mixing is used. 

7.4 Hadronization in BSM models 

In most cases the hadronization of events involving new physics, using the cluster model, proceeds 
in the same way as for Standard Model events. There are however some classes of new physics 
model that require special treatment, in particular: 

Stable strongly interacting particles, if there are strongly interacting particles that are 
stable on the hadronization timescale, these particles will hadronize before they decay. If the 
new particles are in the fundamental representation of colour SU(3) then their hadronization 
proceeds in the same way as for quarks, however if they are in the octet representation the 
situation is more complicated [80]. 

Baryon number violation (BNV) , there are models of new physics in which the conservation 
of baryon number is violated. This typically occurs at a vertex that has the colour tensor e l ^ k 
leading to three quarks, or antiquarks, that are colour connected to each other after the parton 
shower and gluon splitting. 

The Herwig++ hadronization module is designed so that both stable coloured particles and 
baryon number violation are correctly treated as described below. 



w 



w 



w, 



w. 



Ul 



dd 



w 



,1,1 



— COS 

2 

1 • 2 

- sin 

2 



73 



7.4.1 Stable strongly interacting particles 

Currently only the hadronization of objects in the fundamental representation of the SU(3) group 
of the strong force is supported. Provided that the relevant hadrons exist the hadronization of 
these particles is handled in the same way as for quarks. In the future we will extend this to new 
particles in the octet representation as described in Ref. [80]. 

7.4.2 Baryon number violation 

The treatment of QCD radiation and hadronization in models that violate baryon number con- 
servation is described in Refs. [67] and [68] and was implemented in the FORTRAN HERWIG 
program. In events where baryon number is violated there are typically two situations that can 
arise. 

1. The baryon number violating vertices are unconnected, leading to three quarks, or anti- 
quarks, connected to each BNV vertex after the gluon splitting. These (anti)quarks must 
be formed into a cluster, which decays to give a (anti)baryon and a meson, giving the 
expected baryon number violation. In the approach of Refs. [67, 68] this is modelled by 
first combining two of the (anti) quarks into a (anti)diquark, which is in the (anti)-triplet 
representation of colour SU(3). The (anti)quark and (anti)diquark can then be formed into 
a colour singlet cluster, which can be handled by the hadronization module in the normal 
way. 

2. Two baryon number violating vertices are colour connected to each other, leading to two 
quarks connected to one vertex and two antiquarks connected to the second, after gluon 
splitting. In this case two clusters must be formed by pairing one of the quarks with one 
of the antiquarks at random and then pairing up the remaining pair. 

The handling of these colour flows in both the shower and hadronization is fully supported 
although there are currently no models that include baryon number violation implemented. 

7.5 Code structure 

The ClusterHandronizationHandler inherits from the HadronizationHandler of ThePEG and imple- 
ments the cluster hadronization model. The ClusterHandronizationHandler makes use of a number 
of helper classes to implement different parts of the model. The helper classes, in the order they 
are called, are: 

PartonSplitter The PartonSplitter performs the non-perturbative splitting of the gluons in 
quark-antiquark pairs. 

ClusterFinder The ClusterFinder is responsible for taking the partons after the gluon splitting 
and forming them into colour singlet clusters as Cluster particles. 

ColourReconnector It is possible that rather than using the leading N c colour structure of 
the event there is some rearrangement of the colour connections. The option of implementing 
such a model in a class inheriting from the ColourReconnector class is available, although the 
ColourReconnector itself does not implement such a model. 



74 



ClusterFissioner The ClusterFissioner class is responsible for splitting large mass clusters into 
lighter ones as described in Sect. 7.2. 



LightClusterDecayer The LightClusterDecayer decays any clusters for which the decay to two 
hadrons is kinematically impossible into the lightest hadron with the correct flavour together with 
the reshuffling of momentum with neighbouring clusters, which is required to conserve energy 
and momentum, as described at the end of Sect. 7.3. 

ClusterDecayer The ClusterDecayer decays the remaining clusters into pairs of hadrons as 
described in Sect. 7.3. 

In addition to these classes the ClusterDecayer makes use of a HadronSelector to select the 
hadrons produced in the cluster decay 29 . In order to support the different options described 
in Sect. 7.3 the base HadronSelector implements much of the functionality needed to select the 
hadrons in the cluster model but the chooseHadronPair() method, which is used to select the 
hadrons, is virtual and must be implemented in inheriting classes that implement specific variants 
of the cluster model. The FORTRAN HERWIG algorithm is implemented in the Hw64Selector class 
and the Kupco and Herwig++ methods in the HwppSelector class. 

There are a number of switches and parameters that control the hadronization. Here we 
merely give a summary of the most important ones. All the parameters are described in full in 
the Doxygen documentation of the relevant classes. 

The main choice is which variant of the cluster model to use. This can be controlled by 
using either the Hw64Selector for the original model of Ref. [2] or the HwppSelector class for 
the Kupco and Herwig++ variants. The choice of whether to use the Hw64Selector or Hw- 
ppSelector is controlled by setting the HadronSelector interface of the ClusterDecayer and 
LightClusterDecayer classes. In addition, for the HwppSelector the Mode switch controls whether 
the Kupco [Mode=0] or 

Herwig++ [Mode=l], the default, variant is used. The production of specific hadrons by the 
cluster model can be forbidden via the Forbidden interface of the HadronSelector: this option is 
currently only used to forbid the production of the a and k resonances, which are only included 
in the simulation to model low-mass s-wave irn and Kir systems in certain particle decays. 

In addition the mixing angles for the various multiplets can be changed in the HadronSelector 
as can the suppression weights for different SU(3) meson flavour multiplets. 

If the option of using the soft underlying event model [81] is used, as described in Sect. 8.4, 
then the UnderlyingEventHandler needs to be set to the UA5Handler, by default this is set to 
the NULL pointer and the multiple scattering model of the underlying event described in Sect. 8 
used. 

The other main parameters of the cluster model, and their default values, are given in Table 10. 

Finally the ConstituentMass of the gluon and, to a lesser extent the light quarks, which 
can be set in their ParticleData objects, have a major effect on the hadronization since they set 
the scale for the cluster mass distribution. 



3 The LightClusterDecayer also makes use of this class to select the lightest hadron with a given flavour. 
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Parameter 


Default 


Allowed 


Description 




Value 


Range 




HadronSelector 


PwtDquark 


1. 


0-10 


Weight for choosing a down quark 


PwtUquark 


1. 


0-10 


Weight for choosing a up quark 


PwtSquark 


0.68 


0-10 


Weight for choosing a strange quark 


PwtDIquark 


0.49 


0-10 


Weight for choosing a diquark 


SngWt 


0.77 


0-10 


Weight for singlet baryons 


DecWt 


0.62 


0-10 


Weight for decuplet baryons 


LightClusterDecayer 


SingleHadronLimitBottom 


0.16 


0-10 


Bottom cluster to 1 hadron param. 


SingleHadronLimitCharm 


0.0 


0-10 


Charm cluster to 1 hadron param. 


ClusterDecayer 


CIDirLight 


1 


0/1 


Orientation of light cluster decays 


CIDirBottom 


1 


0/1 


Orientation of bottom cluster decays 


CIDirCharm 


1 


0/1 


Orientation of charm clusters 


CISmrLight 


0.78 


0-2 


Smearing of light cluster decays 


CISmrBottom 


0.10 


0-2 


Smearing of bottom cluster decays 


CISmrCharm 


0.26 


0-2 


Smearing of charm cluster decays 


OnShell 





0/1 


Masses of produced hadrons 


ClusterFissioner 


CIMaxLight 


3.15 


0-10 


Max. mass for light clusters (GeV) 


CIMaxBottom 


3.10 


0-10 


Max. mass for bottom clusters (GeV) 


CIMaxCharm 


3.00 


0-10 


Max. mass for bottom clusters (GeV) 


CIPowLight 


1.28 


0-10 


Mass exponent for light clusters 


CIPowBottom 


1.18 


0-10 


Mass exponent for bottom clusters 


CIPowCharm 


1.52 


0-10 


Mass exponent for charm clusters 


PSplitLight 


1.20 


0-10 


Splitting param. for light clusters 


PSplitBottom 


1.00 


0-10 


Splitting param. for bottom clusters 


PSplitCharm 


1.18 


0-10 


splitting param. for charm clusters 


RemnantOption 


1 


0/1 


Treatment of remnant clusters 


SoftClusterFactor 


1 


0.1-10 


Remnant mass param. (GeV) 


ConstituentMasses of light partons (set 


in their ParticleData objects) 


gluon 


0.9 


0-1 


Gluon constituent mass (GeV) 


up 


0.325 


0-m g /2 


Up quark constituent mass (GeV) 


down 


0.325 


0-m s /2 


Down quark constituent mass (GeV) 


strange 


0.5 


m g /2-l 


Strange quark constituent mass (GeV) 



Table 10: Important hadronization parameters. For all parameters other than the light parton 
constituent masses, the limits given are enforced by the interface. For the light partons, the limits 
are not enforced but give a sensible range over which the parameters can be varied. For the gluon, 
the upper limit we give is about the largest value we would consider reasonable, although it is 
not a hard limit. The up and down masses must be less than half the gluon mass, otherwise the 
non-perturbative gluon decays are impossible, and the strange mass must be large enough that 
gluon decays into strange quarks are not possible, to give good agreement with LEP data. 
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8 Underlying Event and Beam Remnants 



The default underlying event model of Herwig++ is currently based on the eikonal model dis- 
cussed in Refs. [18,82-84]. It models the underlying event activity as additional semi-hard and 
soft partonic scatters. In doing so, it allows the description of minimum bias events as well as 
the underlying event in hard scattering processes. The perturbative part of the models provides 
very similar functionality to FORTRAN HERWIG + JIMMY with some minor improvements. The 
non-perturbative part has been newly introduced and is the first available implementation of this 
model. It contains additional developments initiated by the findings in Ref. [85]. 

In this section, we briefly discuss the basics of how to calculate the multiplicities of the semi- 
hard scatterings, before mentioning the details of the soft additional scatters and explaining the 
integration into the full Monte Carlo simulation. For historical reasons, we also briefly mention 
an alternative underlying event model available in Herwig++: the UA5 model [81], even though 
this is ruled out by data and not recommended for serious use. Finally we will describe the code 
structure, which implements these ideas. A more detailed explanation can be found in Ref. [8]. 

8.1 Semi-hard partonic scatters 

The starting point is the observation that the cross section for QCD jet production may exceed 
the total pp or pp cross section already at an intermediate energy range and eventually violates 
unitarity. For example, for QCD jet production with a minimum p t of 2 GeV this already happens 
at ^/s ~ 1 TeV. This pt cutoff should however be large enough to ensure that we can calculate 
the cross section using pQCD. The reason for the rapid increase of the cross section turns out to 
be the strong rise of the proton structure function at small x, since the x values probed decrease 
with increasing centre of mass energy. This proliferation of low x partons may lead to a non- 
negligible probability of having more than one partonic scattering in the same hadronic collision. 
This is not in contradiction with the definition of the standard parton distribution function as the 
inclusive distribution of a parton in a hadron, with all other partonic interactions summed and 
integrated out. It does, however, signal the onset of a regime in which the simple interpretation of 
the pQCD calculation as describing the only partonic scattering must be unitarized by additional 
scatters. 

In principle, predicting the rate of multi-parton scattering processes requires multi-parton 
distribution functions, about which we have almost no experimental information. However, the 
fact that the standard parton distribution functions describe the inclusive distribution gives a 
powerful constraint, which we can use to construct a simple model. The eikonal model used 
in Refs. [18, 82, 83] derives from the assumption that at fixed impact parameter, b, individual 
scatterings are independent and that the distribution of partons in hadrons factorizes with respect 
to the b and x dependence. This implies that the average number of partonic collisions at a given 
b value is 



and cr^ d is the inclusive cross section to produce a pair of partons with p t > p™ m . We model 
the impact parameter dependence of partons in a hadron by the electromagnetic form factor, 





where A(b) is the partonic overlap function of the colliding hadrons, with 




(8.2) 
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resulting in an overlap function for pp and pp collisions of 

A(b^) = -^-^bfK^b), (8.3) 

where /i is the inverse proton radius and K 3 (x) is the modified Bessel function of the third kind. 
We do not fix fi at the value determined from elastic ep scattering, but rather treat it as a free 
parameter, because the spatial parton distribution is assumed to be similar to the distribution 
of charge, but not necessarily identical. 

The assumption that different scatters are uncorrelated leads to the Poissonian distribution 
for the number of scatters, n, at fixed impact parameter, 

V n (b,s)= (n(t Y ))n e ~WM> . (8 .4) 

Using Eq. (8.4) the unitarized cross section can now be written as 

/oo „ 
d 2 b 5> fe (M) = / d 2 b [1- e-^ b ' s »] , (8.5) 
k=i 

which properly takes multiple scatterings into account. The key ingredient for the Monte Carlo 
implementation is then the probability of having n scatterings given there is at least one, inte- 
grated over impact parameter space. This expression reads 

f d 2 b VJb,s) , , 

J d b 22k=i V k{b,s) 

It is worth noting that this distribution, after integration over b, is much broader than Poissonian 
and has a long tail to high multiplicities. 

Equation (8.6) is used as the basis of the multi-parton scattering generator for events in which 
the hard process is identical to the one used in the underlying event, i.e. QCD 2 — ^ 2 scattering. 
For scatterings of more than one type of hard process, the formulae can be easily generalized, 
but in fact for the realistic case in which all other cross sections are small compared to the jet 
cross section, they saturate at a simple form, 



P ' 



"hard 



Jd 2 bV n (b,s), (8.7) 



which allows for a more efficient generation of additional scatterings. It is worth noting that 
the fact that we have 'triggered on' a process with a small cross section leads to a bias in the b 
distribution and hence a higher multiplicity of additional scatters than in the pure QCD 2^2 
scattering case. A slight further modification to the distribution is needed when the small cross 
section process is a subset of the large one, for example QCD 2^2 scattering restricted to the 
high p t region. 

As described so far, the n scatters are completely independent, which is expected to be a good 
approximation in the region in which multiple scattering dominates, i.e. small momentum frac- 
tions. However, some fraction of events come from higher x values and must lead to correlations 
between the scatters at some level. At the very least, the total momentum and flavour must be 
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conserved: the total x value of all partons extracted from a hadron cannot exceed unity and each 
valence parton can only be extracted once. In Herwig++ these correlations are included in the 
simplest possible way, by vetoing any scatters that would take the total extracted energy above 
unity and by only evolving the first scatter back to a valence parton and all the others back to 
a gluon. 



8.2 Soft partonic scatters 

The elastic scattering amplitude, a(h,s), in impact parameter space can be expressed in terms 
of a real eikonal function x(b, s), as 



t (b, S ) = i[e-^ -1] . 



The elastic scattering amplitude, A(s, t), is the Fourier transform of a(b, s) and therefore the total 
PP {pp) cross section as well as the elastic cross section can be obtained from that parameterization 

as, 

a tot (s) = 2 J d 2 b [1 - e"*^ ] , a el {s) = J d 2 b |l - e"*^ f . (8.9) 
The inelastic cross section is obtained as the difference between the two cross sections, 

Cinel = Otot — Cel 



J d 2 b [1 - e" 2 *^ ] . 



(8.10) 



The elastic t-slope parameter at zero momentum transfer is also calculable within this framework 
and yields [86] 

b cl = — [ d 2 b b 2 [l - e"* (b ' s) 1 . (8.11) 

Ctot J 

To reproduce the results from Eq. (8.5), we choose 

X(b )fl ) = i<n(b ja )>. (8.12) 

However we want to introduce additional scatters below the transverse momentum cut-off. 
Therefore, we identify this as the hard part of a universal eikonal function, which then has the 
form, 

Xtot (b,s) = XQCD(b,s) + 

Xsoft 

(b,s), (8.13) 

with the perturbative part 

XQCD (b,s) = ±A(b;n) C d ( S ;pf), (8-14) 

as in Eq. (8.1). 

In the models of Refs. [82-84], the soft eikonal function has the form 

Xsoft(b, s) = ^A soft (b; /i soft ) cr^ft , (8.15) 
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where <r™f t is the purely non-perturbative cross section below pf" n , which is a free parameter 
of the model. That is, we assume that soft scatters are the result of partonic interactions 
that are local in impact parameter. Previous Monte Carlo implementations used the simplest 
assumption about the partonic overlap function probed by the soft scatters, v4 so f t (b) = A(h), 
i.e. an identical distribution to the one probed by semi-hard scatters. In Ref. [85] it was shown 
that measurements on the elastic t-slope confine the allowed parameter space of such models 
vastly. The remaining parameter space seems to be in contradiction with constraints obtained 
from measurements of the effective cross section in double parton scattering events [87, 88]. 
We therefore introduced the option of relaxing the condition of identical overlap distributions in 
Herwig++. The option twoComp of the MPIHandler enables the dynamical determination of the 
soft overlap distribution, A so f t (b). In this case, which is the default setting, we use Eq. (8.3) but 
allow an independent radius parameter for the soft overlap function. The parameter /i so f t is then 
dynamically fixed by the requirement of a correct description of the elastic t-slope from Eq. (8.11) 
at the current centre-of-mass energy. At the same time we fix the second free parameter in the soft 
sector, <7™f t , by choosing it such that the total cross section, evaluated with the parametrization 
from Ref. [89] is correctly described. First measurements of the total cross section may deviate 
from the prediction in Ref. [89] and therefore the parameter MeasuredTotalXSec can be used 
to set the total cross section at the current centre-of-mass energy explicitly. 

With the full eikonal from Eq. (8.13), we can construct our model for additional semi-hard 
and soft scatters, by imposing the additional assumptions, 

• The probability distributions of semi-hard and soft scatters are independent 

• Soft scatters are uncorrelated and therefore obey Poissonian statistics like the semi-hard 
scatters 

The probability 7\ n (b, s), for having exactly h semi-hard and n soft scatters at impact parameter 
b and centre-of-mass energy s is then given by, 



From Eq. (8.16) we can now deduce the cross section for having exactly h semi-hard and n soft 
scatters as, 



The cross section for an inelastic collision (either semi-hard or soft), is obtained by summing over 
the appropriate multiplicities and yields 








h+n>l 





The inelastic cross section for at least one semi-hard scattering is 




/i>l,n>0 




(8.19) 
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With the cross sections from Eqs. (8.17) and (8.18) we can construct the basis of our multiple 
soft and semi-hard scattering model, the probability, Ph, n , of having exactly h semi-hard and n 
soft scatters in an inelastic event (h + n > 1). It is given by 

p , s °h,n( s ) Jd 2 b 7\ w (b,s) 

Phn(s) — —TT = r F 7T — rr , fl + U > 1 , 8.20 

fTinel(s) / d 2 b [1 - e-Xt°tal(b, S ) ] ' - ' v ) 

which is analogous to Eq. (8.6) for the case of solely semi-hard additional scatterings. Equa- 
tion (8.20) defines a matrix of probabilities for individual multiplicities. This matrix is evaluated 
at the beginning of each run and the corresponding multiplicities are drawn for each event from 
this matrix according to their probability. 

Equation (8.20) leads to very inefficient generation of additional scatters in cases where a rare 
hard scattering, with cross section a iauTC , takes place. Equation (8.7) has been deduced for this 
case, by exploiting the independence of different scatters. The presence of soft scatters does not 
alter that result as our assumption is that the soft scatters are independent from each other and 
from the other scatterings. Hence, the probability for h hard scatters (from which one is distinct, 
i.e. h = m + 1) and n soft scatters is given by 

. Jd 2 hV mn (h,s) W*)^) 1 e -m^e 
Ph - " '/d'b A{b)a^ c (8 ' 21) 



J d 2 bP ra , n (b,^(6) (8.22) 
d 2 b V h Jb, s) . (8.23) 



h 



'hard 



The probability for m semi-hard (pt > pf 1111 ) and n soft additional scatters is multiplied in 
Eq. (8.21) with the probability for exactly one scattering with an inclusive cross section of o" rarc . 
The denominator is the inclusive cross section for this distinct scattering, i.e. summed over all 
multiplicities for additional semi-hard and soft scatters. By approximating the exponential with 
unity and exploiting the normalization of A(b) (J d 2 b A(b) = 1), we finally deduce Eq. (8.23). 



8.2.1 Monte Carlo implementation 

In this section, we describe in detail how the additional soft scatterings are implemented. The 
corresponding description for the semi-hard part of the underlying event is given in Ref. [8]. 

At large centre-of-mass energies, s, and small scale of interactions, Q 2 , the parton densities 
suggest a proliferation of small-x gluons. That is the reason why we chose to model the soft 
scattering contributing to the underlying event as elastic collisions between soft gluons. We start 
the generation of these soft scatters after all perturbative evolution has terminated, since p™ m is 
typically at the order of the parton shower cut-off scale. The non-perturbative remnant decays, 
that are described in detail in Ref. [8] , produce diquarks from which the soft gluons are radiated 
and scatter off each other. Such a scattering is depicted in Fig. 5. 

All soft gluons carry colour charge and have an effective gluon mass 30 , m g = 0.75 GeV, in 
correspondence to the effective gluon mass that is used during parton showers and hadronization. 

30 This is currently hard coded, but could be linked to the gluon mass used by the hadronization model in future 
versions. 
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Figure 5: Soft gluon collision in a diquark scattering. The diquarks are in a anti-triplet state 
and remain unchanged with respect to their colour state. 



As the simplest solution, we have chosen to sever the colour connections to the diquarks so that 
the two outgoing gluons from each soft scattering are colour connected to each other, similar to 
the colourDisrupt=l option for the semi-hard scattering. This is motivated by a Pomeron-like 
structure for these soft forward interactions. 

The scattering of the soft gluons can be described by the variables x\,X2,Pt, 4>- %i and x 2 
are the longitudinal momentum fractions of the two incoming gluons (<7i,<72)j so that their 4 
momenta in the lab frame are 




x 



p t and are the transverse momentum and azimuthal angle, in the CM frame, of the outgoing 
gluons (#3, g^) respectively. Their four-momenta can therefore be parameterized as, 

Pg 3A = (y\Jp 2 t+P 2 z + m l , ±p t cos0, ±p t sm<f), ±p^j , (8.25) 
where the longitudinal momentum, p z , is fixed by total energy-momentum conservation, 

Pz = 4 Pt-m 9 - (8-26) 

The kinematics of the soft processes are fixed by choosing values for the four parameters. x\ 
and x 2 are sampled from a f(x) = 1 /^-distribution in the range [x min ,x max ]. x min is a cut-off to 
avoid numerical instabilities 31 . x max corresponds to the maximum available energy that is left 
in the diquarks. The azimuthal angle is sampled from a uniform distribution, G (0,27r). The 
transverse momentum is the last remaining degree of freedom. By construction the transverse 
momentum distribution must not exceed p™ m , but the functional form of it is not predetermined. 
We use a Gaussian distribution, 

( ^ = Ae~^ ) (8.27) 
dpf 

to parameterize it. To fix the free parameters A and (3, we impose the following constraints: 



31 



At present this is hard coded as x m i n = (2p™ ln ) js 
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Figure 6: Transverse momentum distribution of additional scatters 



The resulting soft cross section has to match the total soft cross section, which has been 
fixed to describe <7 to t and b e \ with Eqs. (8.9) and (8.11), 

H (T inC I 

M = <ft- (8-28) 

The transverse momentum distribution of semi-hard and soft scatterings should be contin- 
uous at the matching scale, 



J inc 

#(s;p™ in ) :- hard 



rl rr lnc 
ucr soft 



dp} 



where we introduced H as shorthand for the hard inclusive jet cross section at p t = pi 

These conditions are fulfilled by the parameterization, 

da, 
dpj 

where the slope, /3, must satisfy, 



(8.29) 



^ = H(s; pf n ) e -^- p ? in ) , (8.30) 



,0pT in2 - 1 rr inc 

soft (8.3r 



P H(s;p?™) 

Figure 6 shows the transverse momentum spectrum for two different cut-off values. The slope, 
(3, is chosen such that both curves correspond to the same integrated cross section. 

After the kinematics have been generated in the CM frame, we boost back to the lab frame 
and reshuffle the diquark momenta such that they remain on their original mass shell. Now we 
can determine the available energy for the next soft interaction and iterate the process until the 
requested multiplicity has been reached or all available energy of the diquarks has been used. 
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8.3 Connection to different simulation phases 

The model introduced so far is entirely formulated at the parton level. However, an event 
generator aims for a full description of the event at the level of hadrons. This implies that the 
implementation of multi-parton scattering must be properly connected to the parton shower and 
hadronization models, a few details of which we discuss in the following. 

8.3.1 Parton showers and hard matrix elements 

After generating the hard process and invoking parton showers on its coloured particles, the 
number of additional scatters is calculated according to Eq. (8.20) or Eq. (8.23) respectively. 
After the initial-state shower has terminated, the incoming partons are extracted out of the 
beam particles in the usual way. 

The requested additional scatters are then generated using a similar but completely inde- 
pendent infrastructure from the one of the hard process. Dedicated hard matrix elements with 
hand-coded formulae summed over parton spins are used for greater speed, as mentioned in 
Sect. 3.1. This also has the advantage that specific cuts, different to those used for the main 
hard process, can be specified. 

For each additional scattering, parton showers evolve the produced particles down to the 
hadronic scale. The backward evolution of additional scatters is forced to terminate on a gluon. 
After termination, these gluons are extracted out of the beam particles. If this process leads to 
a violation of four-momentum conservation, the scattering cannot be established. It is therefore 
regenerated until the desired multiplicity has been reached. If a requested scattering can never 
be generated without leading to violation of momentum conservation, the program eventually 
gives up, reducing the multiplicity of scatters. 

8.3.2 Minimum bias process 

Starting from version 2.3, Herwig++ simulates minimum bias collisions as events in which there 
is effectively no hard process. However, to maintain a uniform structure with the simulation 
of standard hard process events, we have implemented a matrix element class, MEMinBias that 
generates a 'hard' process with as minimal an effect as possible. It extracts only light (anti) quarks 
(d, u, d or u) from the hadrons and allows them to 'scatter' through colourless exchange at zero 
transverse momentum, with matrix element set to unity, so that their longitudinal momentum 
is determined only by their parton distribution functions. To give a predominantly valence-like 
distribution, a cut on their longitudinal momentum fraction x > 10~ 2 is recommended, as shown 
in Sect. B.2. Note that because the matrix element is set to unity, the cross section that is printed 
to the output file at the end of the run is meaningless. 

8.3.3 Hadronization 

The underlying event and beam remnant treatment are closely connected because the generation 
of additional scatters requires the extraction of several partons out of the proton. As described 
before, all additional partons are extracted from the incoming beam particles. This is different 
from the procedure that was used in FORTRAN JIMMY, where the successive partons were always 
extracted from the previous beam remnant, a difference in the structure of the event record that 
should not lead to significant differences in physical distributions. 
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Name 


Description 


Default 


Nl 
N2 
N3 


a in n = a (s/GeV 2 ) & + c 

/) in T) — n ( q i ' C^rC*\l'^'\ —1— r> 
U ill / i — ti \o / vi i v J 

c in n = a (s/GeV 2 ) 6 + c 


9.110 

yj. J. ± o 

-9.500 


K 1 

K2 


a m = am Is/ Lre V ) + o 
b in 1/Jfe = a In (s/GeV 2 ) +6 


n non 
O.Ozy 

-0.104 


Ml 

M2 


a m — mi — m 2 — aje 
6 in (M — mi — m 2 — a)e~ bM 


0.4 GeV 
2.0 GeV" 1 


PI 
P2 
P3 


Pt slope for d, u 
Pt slope for s,c 
p t slope for qq 


5.2 GeV" 1 
3.0 GeV" 1 
5.2 GeV" 1 



Table 11: Parameters of the soft underlying event event model. 

The cluster hadronization described in the previous section can only act on (anti) quarks or 
(anti)diquarks. However, naively extracting several partons from a hadron will not leave a flavour 
configuration that is amenable to such a description in general. Therefore, the strategy we use, 
as already mentioned, is to terminate the backward evolution of the hard process on a valence 
parton of the beam hadron and additional scatterings on gluons, giving a structure that can be 
easily iterated for an arbitrary number of scatters. This structure is essentially the same as in 
FORTRAN JIMMY. 

8.4 UA5 parametrization 

While the new multiple interaction model provides a better description of the underlying event 
and is recommended for all realistic physics studies, Herwig++ still contains the original soft 
model of the underlying event used before version 2.1. 

This model is based on the minimum-bias event generator of the UA5 Collaboration [81], 
which starts from a parameterization of the pp inelastic charged multiplicity distribution as a 
negative binomial distribution. In Herwig++ the relevant parameters are made available to the 
user for modification, the default values being the UA5 ones as used in the FORTRAN version of 
the program. These parameters are given in Table 11. 

The first three parameters control the mean charged multiplicity n at cm. energy y/s as 
indicated. The next two specify the parameter k in the negative binomial charged multiplicity 
distribution, 

r(n + k) (n/kr 
{ ' n\V{k) (l + n/k) n + k ' 

The parameters Ml and M2 describe the distribution of cluster masses M in the soft collision. 
These soft clusters are generated using a flat rapidity distribution with Gaussian shoulders. The 
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transverse momentum distribution of soft clusters has the form 

P(p t ) oc p t exp \ -by/j% + M 2 ) , 

where the slope parameter b depends as indicated on the flavour of the quark or diquark pair 
created when the cluster was produced. As an option, for underlying events, the value of y/s 
used to choose the multiplicity n may be increased by a factor EnhanceCM to allow for an 
enhanced underlying activity in hard events. The actual charged multiplicity is taken to be n 
plus the sum of the moduli of the charges of the colliding hadrons or clusters. 

8.5 Code structure 

In addition to being the main class responsible for the administration of the shower, the 
ShowerHandler, described in Sect. 6.10, is also responsible for the generation of the additional 
semi-hard scattering processes. It has a reference to the MPIHandler set in the input files, which 
is used to actually create the additional scattering processes. It invokes the parton shower on all 
the available scatters and connects them properly to the incoming beam particles. This includes 
potential re-extraction of the incoming parton if it is changed as a result of initial-state radiation. 
It modifies the RemnantParticles that were initially created by the Parton Extractor. A number 
of classes are used by the ShowerHandler to generate the additional scattering processes. Soft 
additional scatters are generated in the HwRemDecayer class. 

MPIHandler The MPIHandler administers the calculation of the underlying event activity. It 
uses MPISampler to sample the phase space of the processes that are connected to it. Using that 
cross section the probabilities for the individual multiplicities of additional scatters are calculated 
during initialization. The method MPIHandler: :multiplicity() samples a number of extra scatters 
from that pretabulated probability distribution. The method MPIHandler::generate() creates one 
subprocess according to the phase space and returns it. 

MPISampler The MPISampler performs the phase-space sampling for the additional scatter- 
ings. It inherits from SamplerBase and implements the Auto Compensating Divide-and-Conquer 
phase space generator, ACDCGen. 

HwRemDecayer The HwRemDecayer is responsible for decaying the RemnantParticles to some- 
thing that can be processed by the cluster hadronization, i.e. (anti)quarks or (anti)diquarks. This 
includes the forced splittings to valence quarks and gluons respectively. Also the colour connec- 
tions between the additional scatters and the remnants are set here. If additional soft partonic 
interactions, i.e. the non-perturbative part of the underlying event, are enabled, they are gener- 
ated inside this class after the remnants have been decayed to the (anti)diquarks. 

MPIPDF The MPIPDF class is used to modify the PDF's used for the initial state shower of ad- 
ditional scatters. All sorts of rescaling are possible but currently the mode that is used is the one 
where the valence part of the PDF is removed. The objects are instantiated inside ShowerHandler 
and set to the default PDF's using void ThePEG: : CascadeHandler : :resetPDFs(. . .) 

The most important interfaces to set parameters for the above mentioned classes are described 
here. An exhaustive description of all interfaces is provided by our Doxygen documentation. 
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MPIHandler 



SubProcessHandlers: Vector of references to SubProcessHandler objects. The first ele- 
ment is reserved for the underlying event process. Additional references can be set to simulate 
additional hard processes in a single collision. See Sect. B.3 for details of how to use this func- 
tionality. 

Cuts: Vector of references to Cuts objects. The first element is used to impose the minimal pt 
of the additional scatters, p™ m . This is one of the two main parameters of the model. The current 
default, obtained from a fit to Tevatron data is 4.0 GeV. See Ref. [19] for details. Additional 
cuts object may be defined for additional hard processes that should be simulated in the same 
event. 

additionalMultiplicities: Vector of integer values to specify the multiplicity of additional 
hard scattering processes in a single collision. See Sect. B.3 for an example. 

InvRadius: The inverse beam particle radius squared, /i 2 . The current default is 1.5 GeV 2 , 
obtained from the above mentioned fit. 

IdenticalToUE: An integer parameter specifying which element of the list of 
SubProcessHandlers in SubProcessHandlers is identical to the underlying event process. Zero 
means the the conventional hard subprocess is QCD jet production. -1 means that no process 
is identical. Any number > means that one of the additional hard scatterings is QCD jet 
production, where the exact number specifies the position in the vector. The default is -1, which 
is appropriate as long as no QCD jet production is simulated. 

colourDisrupt: Real number in the range [0, 1], which gives the probability for an ad- 
ditional semi-hard scattering to be disconnected from other subprocesses as far as the colour 
connections are concerned. The current default is 0. 

softlnt: Switch to turn the inclusion of non-perturbative scatters to the underlying event 
model on (Yes) or off (No). The current default is Yes. 

twoComp: Switch to toggle between an independent overlap function for soft additional 
scatters (Yes) and identical ones A so f t (b) = A(h) (No). If the two-component model is used, the 
parameters of the soft sector are automatically choosen to describe the total cross section as well 
as the elastic t-slope correctly. 

DLmode: Integer number e {1, 2, 3} to choose between three different parametrizations of 
the total cross section as a function of the centre-of-mass energy: 

1. Parametrization of Ref. [89]. 

2. Parametrization of Ref. [89] but with rescaled normalization to match the central value of 
the measurement [90] by CDF. Default 
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3. Parameterization of Ref. [91]. 

MeasuredTotalXSec: Parameter to set the total cross section (in mb) explicitly. If this 
parameter is used, it will overwrite the parametrization selected with the previous switch. This is 
intended for first data on the total cross section and should be used instead of the parametrization, 
which may deviate substantially. 

ShowerHandler 

MPIHandler: Reference to the MPIHandler. To switch multiple parton interactions off, 
this reference has to be set to NULL. 



Since it is not a recommended option, we do not go into as much detail, but for complete- 
ness, we briefly mention the structure of the UA5 code. The UA5 model is implemented in 
the UA5Handler class, which is intended to be used as an UnderlyingEventHandler in the 
ClusterHadronizationHandler. The main interfaces of the UA5Handler are the parameters named 
in Table 11, described in Sect. 8.4. 
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9 Hadron Decays 



Herwig++ uses a sophisticated model of hadronic decays, as described in Refs. [25,92]. The 
simulation of decays in Herwig++ is designed to have the following properties: 

• a unified treatment of the decays of both the fundamental particles and the unstable 
hadrons, this is of particular importance for particles like the r lepton, which, while a 
fundamental particle, is more akin to the unstable hadrons in the way it decays; 

• up-to-date particle properties, i. e. masses, widths, lifetimes, decay modes and branching ra- 
tios together with a new database to store these properties to make updating the properties 
easier and the choices made in deriving them clearer; 

• full treatment of spin correlation effects using the algorithm of Refs. [27-30] for the decay 
of all unstable particles, it is important that the same algorithm is used consistently in 
all stages of the program so that correlations between the different stages can be correctly 
included; 

• a sophisticated treatment of off-shell effects for both the unstable hadrons and fundamental 
particles; 

• a large range of matrix elements for hadron and tau decays including both general ma- 
trix elements based on the spin structures of the decays and specific matrix elements for 
important decay modes; 

• the accurate simulation of QED radiation in the particle decays using the Yennie-Frautschi- 
Suura (YFS) formalism. 

In this section we describe the simulation of hadron and tau decays in Herwig++. We start 
by discussing the physical properties of the hadrons used in the simulation and how they are 
determined. In ThePEG framework these physical properties are stored using the ParticleData 
class, which has one instance for each particle used in the simulation. In turn the properties 
of a given decay mode are stored using the DecayMode class, which contains both the particles 
involved in the decay and a reference to a Decayer object that can be used to generate the 
kinematics of the decay products. The DecayHandler class then uses these DecayMode objects to 
select a decay of a given particle, according to the probabilities given by the branching ratios for 
the different decay modes, and then generates the kinematics using the relevant Decayer specified 
by the DecayMode. 

Following a brief discussion of the treatment of off-shell effects we therefore discuss the differ- 
ent Decayer classes available in Herwig++ for the decay of tau leptons, strong and electromagnetic 
hadron decays and then hadron decays. This is followed by a discussion of the code structure. 

9.1 Particle properties 

The information in the Particle Data Group's (PDG) compilation [49] of experimental data is 
sufficient in many cases to determine the properties of the hadrons used in Herwig++. However, 
there are some particles for which the data are incomplete or too inaccurate to be used. Equally, 
there are a number of particles that are necessary for the simulation but have not been observed, 
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particularly excited bottom and charm hadrons, which should perhaps be regarded as part of the 
hadronization model affecting the momentum spectrum of lighter states, rather than as physical 
states. A large number of choices therefore have to be made in constructing the particle data 
tables used in the event generator based on the data in Ref. [49]. 

In the past the data were stored as either a text file or the contents of a FORTRAN COMMON 
block. However, due to the relatively large amount of data that needs to be stored we decided 
to adopt a database approach based on the MySQL database system. The event generation still 
uses text files to read in the particle properties but these files are now automatically generated 
from the database. This provides us with a range of benefits: the data can now be edited using 
a web interface; additional information describing how the particle properties were determined 
is stored in the database both improving the long-term maintenance and allowing the user to 
understand the uncertainties and assumptions involved. 

An example of the output from the database for the properties of the to meson is shown 
in Figure 7. This includes the basic properties of the particle together with an explanation of 
how they were derived. In addition there is a star rating between one and five, which gives an 
indication of how reliable the properties of the particle and the modelling of individual decay 
modes are. 

In general we used the following philosophy to determine the particle properties used in 
Herwig++: 

• The properties of the light mesons in the lowest lying multiplets were taken from Ref. [49] . 
In some cases we used either lepton universality or the phase-space factors from our Decayers 
to average the branching ratios for poorly measured modes. 

• Where possible the properties of the excited light mesons were taken from Ref. [49] together 
with some additional interpretation of the data. Except for the l 3 Di multiplet, which is 
missing a 0-like member, the mesons needed to fill the l 1 So, l 3 Si, l^Pi, l 3 Po ; l 3 Pi, 1 3 P2, 
1^2, l 3 Di, 1 3 D 3 , 2%, 1% and 2 3 Si SU(3) multiplets are included together with the K 
mesons from the 1 3 D 2 multiplet. 

• The properties of the D u ^ s mesons were taken from Ref. [49] together with the addition of 
some high multiplicity modes to ensure that the branching ratios sum to one. 

• The branching ratios and properties for B u ^ s mesons were taken from the data tables 
of EvtGEN [93], which have been extensively tuned to 5-factory data. This means that 
partonic decay modes are used to model many of the inclusive B decay modes. 

• The mass of the B c meson is taken from Ref. [49] . The branching ratios were taken from 
the theoretical calculations of Ref. [94] together with some partonic modes to ensure that 
the branching ratios sum to one. 

• The properties and decay modes of the charmonium resonances were taken from Ref. [49] 
where possible together with the use of partonic decays, to ggg, gg or qq, to model the 
unobserved inclusive modes. For some of the particles, in particular the h c and rj c (2S), 
the results of Ref. [95] were used and the rj c (2S) branching ratios were taken from the 
theoretical calculation of Ref. [96] . 
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The uj is the lightest isospin singlet from the l 3 5i multiplet. The modes and properties are taken from Ref. [49] 
with the lepton modes averaged. The modes with photons that can be produced by QED radiation are included 
in the mode without the radiation. The uj is allowed to be off-shell by ten times the width. The uj has mass 
782.65 McV and is unstable. The uj has spin 1, charge and is colour neutral. The uj is a meson and is from 
the l 3 Si multiplet. The uj has width 8.49 McV. The lower limit on the mass of the particle is 84.9 MeV and the 
upper limit is 84.9 MeV. These are the deviations from the on-shell value. The branching ratios are fixed. The 
PDG code is 223. The mass generator is omegamass for the uj. The width generator is omegawidth for the uj. 
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Figure 7: An example of the particle properties in Herwig++, in this case for the u meson. The 
properties of the particle including the mass, width, decay modes and branching ratios are given together 
with comments on how properties were determined. In the full web version links are included to the 
descriptions of the objects responsible for generating the kinematics for the various decay modes. 

• The properties and decay modes of the bottomonium resonances were taken from Ref. [49] 
where possible. In addition we have added a large number of states that are expected 
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to have small widths, i.e. the mass is expected to be below the BB threshold, using the 
theoretical calculations of Refs. [97-101] for many of the properties. 

The properties of the excited D and D s mesons were taken from Ref. [49] including recent 
results for the D[ and Dq states. The widths of the D s i and D s2 mesons were from the 
theoretical calculations of Ref. [102] and Ref. [103], respectively. For many of the mesons 
we were forced to assume that the observed modes saturated the total width in order to 
obtain the branching ratios using the results in Ref. [49]. 

The properties of the excited B u ^ s mesons are uncertain. The B* u d s have been observed 
and there is evidence in Ref. [49] from LEP for further excited states, however it was unclear 
which states have been observed. There have been recent claims for the observation of the 
B x , B* and B* 2 states by CDF and DO [104, 105] and the B al by CDF. The situation is still 
unclear, the masses measured by the two experiments disagree for the Bi, B 2 states and DO 
do not observe the B s \ state. We have chosen to use the DO results for the B system and 
the CDF results for the B s system for the observed states and have taken the properties of 
the remaining unobserved states from Ref. [103]. 

The masses of the excited B c mesons, which have not been observed, are taken from the 
lattice results in Ref. [106] , which agree with potential model calculations. The widths and 
branching ratios were taken from the theoretical calculation of Ref. [107]. 

The properties of the light baryons were taken from Ref. [49] where possible. In general 
we have included all the light baryons from the first (56, 0q ) octet and decuplet multiplets. 
We now include the light baryons from the next ^ (56, 0^~), ^ (70, 1^ ), and | (70, 1^ ) 
multiplets, although in some cases we have used higher H resonances whose properties are 
better determined rather than those given in Ref. [49]. In addition the singlet A(1405) 
and A(1520) are also included. By default the | (70, lj~) multiplet and A(1520) are not 
produced in the hadronization stage in order to improved the agreement with LEP data. 

The properties of the weakly decaying charm baryons were taken from [49] together with a 
number of decay modes with theoretical calculated branching ratios from [108] and partonic 
decay modes in order to saturate the total width. 

The experimental data on the weakly decaying bottom baryons is limited. Where possible 
this data, taken from Ref. [49], was used together with a number of theoretical calcula- 
tions [109-114] for the branching ratios to exclusive modes. The masses were calculated 
using the equivalent splitting in the charm system and the A& mass where they have not 
been measured. In addition to the exclusive modes a number of partonic modes are included 
to model the unobserved exclusive decays. 

The properties of the strongly and radiatively decaying charm baryons, i.e E c , S' c , and 
excited A c and S c , are taken from Ref. [49] together with some results from Ref. [115] for 
branching ratios and widths where the experimental data is insufficient. 

The baryons containing a single charm quark from the multiplets containing the A (1405) 
and A(1520) have been observed and are included with the properties taken from Ref. [49] 
where possible and Ref. [115] for some widths. 
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• The same set of excited baryons containing a bottom quark are included as in the charm 
system, despite none of these particles having been observed. The masses are calculated 
using the equivalent splitting in the charm system and the A& mass and the branching ratios 
are assumed to be the same as for the corresponding charm decay. The widths are taken 
from [115]. 

• Given that baryons containing more than one heavy quark can not be produced in the 
cluster hadronization model none of these states, or pentaquarks, are currently included in 
the particle properties. 

All the particle properties used in Herwig++ can be accessed via the online interface to our 
database of particle properties at 

http : / /www . ippp . dur . ac . uk/~richardn/particles/ 



9.2 Line shapes 

In general, if we wish to include the off-shell effects for an outgoing external particle in a hard 
production or decay process we need to include the following factor in the calculation of the 
matrix element 

If mT(m) 
W oS = - dm 2 — K > 2r2 , r , 9.1 

where M is the physical mass of the particle, m is the off-shell mass and T(m) is the running width 
evaluated at scale m. In practice other effects can be included to improve this simple formula, for 
example we include the Flatte lineshape [116] for the ao(980) and /o(980) mesons. In Herwig++ 
we calculate the running width of the particle based on its decay modes. The Decayer responsible 
for each decay mode specifies the form of the running partial width, F^m), for the decay mode 
either in a closed analytic form for two-body decays or as a WidthCalculatorBase object, which 
is capable of calculating the partial width numerically and is used to construct an interpolation 
table. The running width for a given particle is then the sum of the partial widths 



r(m) = 5^(771). (9.2) 



This gives both a sophisticated model of the running width based on the decay modes and allows 
us to use the partial widths to normalize the weights for the phase-space integration of the decays 
to improve efficiency close to the kinematic threshold for the decay. 

In some cases, where the partial width varies significantly over the mass range allowed in the 
simulation, we can choose to use a variable branching ratio 

BRj(m) = ^ (9 .3) 
1 (m) 

both to prevent the production of kinematically unavailable modes and to improve the physics 
of the simulation. The classic examples are the decays of the /o and a$ scalar mesons, which lie 
close to the KK threshold. This means that, depending on their mass, they decay to either -kit 
or r/7r respectively below the threshold or with a significant KK branching fraction above the 
KK threshold. 
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The weight in Eq. (9.1) is automatically included for all the Decayers inheriting from the 
Decay I integrator class, which is the case for vast majority of the Herwig++ Decayers. The 
GenericWidthGenerator calculates the running widths using information from the Herwig++ De- 
cayers inheriting from the Decay I integrator class. For decayers inheriting from the 
BaryonlMesonDecayerBase the running width is calculated using the BaryonWidthGenerator class. 
GenericMassGenerator is responsible for calculating the weight in Eq. (9.1) or generating a mass 
according to this distribution. 

9.3 Tau decays 

The simulation of r lepton decays in Herwig++ is described in detail in Ref. [25], together with 
a detailed comparison between the results of Herwig++ and TAUOLA [117, 118]. Here we simply 
describe the basic formalism for the decays of the tau and the different models available for the 
different decays, together with the structure of the code. 

The matrix element for the decay of the r lepton can be written as 

M = -^L ft J», = u(p Vr ) 7m (1 - 75) u(p T ), (9-4) 

where p T is the momentum of the r and p Vr is the momentum of the neutrino produced in the 
decay. The information on the decay products of the virtual W boson is contained in the hadronic 
current, J M . This factorization allows us to implement the leptonic current L M for the decaying 
tau and the hadronic current separately and then combine them to calculate the r decay matrix 
element. 

In Herwig++ this factorization is used to have a TauDecayer class, which implements the cal- 
culation of the leptonic current for the r decay and uses a class inheriting from WeakDecayCurrent 
to calculate the hadronic current for a given decay mode. This factorization allows us to reuse 
the hadronic currents in other applications, for example in weak meson decay using the naive 
factorization approximation or in the decay of the lightest chargino to the lightest neutralino 
in Anomaly Mediated SUSY Breaking (AMSB) models where there is a small mass difference 
between the neutralino and chargino. 

9.3.1 Hadronic currents 

We have implemented a number of hadronic currents, which are mainly used for the simulation of 
r decays. These are all based on the WeakDecayCurrent class. In this section we list the available 
hadronic currents together with a brief description, a more detailed description can be found in 
either Ref. [25] or the Doxygen documentation. 

ScalarMesonCurrent The simplest hadronic current is that for the production of a pseu- 
doscalar meson, e.g. the current for the production of in the decay of the tau. The hadronic 
current can be written as 

J" = f P P% (9.5) 

where p P is the momentum of the pseudoscalar meson and fp is the pseudoscalar meson decay 
constant. 
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VectorMesonCurrent The current for the production of a vector meson is given by 

r = V2g v e^, (9.6) 

where e*y is the polarization vector for the outgoing meson and gy is the decay constant of the 
vector meson. 

LeptonNeutrinoCurrent The current for weak decay to a lepton and the associated anti- 
neutrino is given by 

J" = 12(^)7^(1-75)^), (9.7) 
where p„ is the momentum of the anti-neutrino and pt is the momentum of the charged lepton. 

TwoMesonRhoKStarCurrent The weak current for production of two mesons via the p or 
K* resonances has the form 

J" = (pi - pa), (sT - ^f) E ^BW k (q 2 ), (9-8) 

V 1 J l^k a k k 

where pi^ are the momenta of the outgoing mesons, q = pi + p 2 , BW k {q 2 ) is the Breit-Wigner 
distribution for the intermediate vector meson k and a.f. is the weight for the resonance, which 
can be complex. The Breit-Wigner terms are summed over the p or K* resonances that can 
contribute to a given decay mode. 

The models of either Kiihn and Santamaria [119], which uses a Breit-Wigner distribution 
with a p-wave running width, or Gounaris and Sakurai [120] are supported for the shape of the 
Breit-Wigner distribution. 

KPiCurrent Unlike the 7r + 7r° decay of the tau the Ktt decay mode can occur via either in- 
termediate scalar or vector mesons. We therefore include a model for the current for the Ktt 
decay mode including the contribution of both vector and scalar resonances based on the model 
of Ref. [121]. The current is given by 

J^ = Cy(p 1 -p 2 ) u ^J—Y,^BW k (q 2 ) (g^-^) +Cs q»=J—J2^ BW ^ 2 ^ ( 9 - 9 ) 

l^k a k £ V M k J Z^k Pk ^ 

where p\^ are the momenta of the outgoing mesons, q = pi +P2, BW*,.(g 2 ) is the Breit-Wigner 
distribution for the intermediate mesons and a k is the weight for the resonance. The sum over 
the resonances is over the vector K* states in the first, vector, part of the current and the excited 
scalar K* resonances in the second, scalar, part of the current. By default the vector part of the 
current includes the i^*(892) and K*(1410) states and the scalar part of the current includes the 
i^g(1430) together with the option of including the k(800) to model any low-mass enhancement 
in the mass of the Kn system, although additional resonances can be included if necessary. 
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ThreeMesonCurrentBase In order to simplify the implementation of a number of standard 
currents for the production of three pseudoscalar mesons we define the current in terms of several 
form factors. The current is defined to be [117] 

^ = (V - ^if) [Fi( P2 - p 3 T + F 2 {p z - PiY + FsiPi ~ P 2 )1 (9.10) 

+<f F 4 + iFse^pfpPpl 

where pi,2,3 are the momenta of the mesons in the order given below and are the form 

factors. We use this approach for a number of three-meson modes that occur in r decays: 
it-it-it + ; TrW"; K~n~K + ; K°7r~K°; K--k K°; vrVif"; K^tt~tt + ; tt~K°tt°; tt~tt°7]; K^Kf, 
K\ts~ K\\ Kg7r~Kl. The current is implemented in terms of these form factors in a base class so 
that any model for these currents can be implemented by inheriting from this class and specifying 
the form factors. 

We currently implement three models for these decays, the 
ThreeMesonDefaultCurrent model of Refs. [117,119,122], which treats all the decay modes, the 
ThreePionCLEOCurrent model of CLEO [123] for the three-pion modes and the 
KaonThreeMesonCurrent model of Ref. [124] for the kaon modes. 

ThreeMesonDefaultCurrent This is the implementation of the model of Refs. [117, 119, 122] , 
which uses the form of Ref. [119] for the a± width. The form factors for the different modes are 
given in Refs. [117,122]. 

ThreePionCLEOCurrent This is the implementation of the model of Ref. [123] for the weak 
current for three pions. This model includes p mesons in both the s- and p-wave, the scalar a 
resonance, the tensor f2 resonance and scalar /o(1370). The form factors for the 7r°7r°7r~ mode 
are given in Ref. [123] and the others can be obtained by isospin rotation. 

KaonThreeMesonCurrent Like the model of Ref. [122] the model of Ref. [124] is designed 
to reproduce the correct chiral limit for tau decays to three mesons. However, this model makes 
a different choice of the resonances to use away from this limit for the decays involving at least 
one kaon and in the treatment of the K\ resonances. The form factors for the different modes 
are given in Ref. [124]. 

TwoPionPhotonCurrent The branching ratio for the decay r~ — > lot\~v t is 1.95% [49]. The 
majority of this decay is modelled as an intermediate state in the four-pion current described 
below. However there is an 8.90% [49] branching ratio of the io into 7r°7, which must also be 
modelled. We do this using a current for 7r ± 7r°7 via an intermediate u. The hadronic current for 
this mode, together with the masses, widths and other parameters, are taken from Ref. [117]. 

FourPionNovosibirskCurrent We use the model of Ref. [125] 32 to model the decay of 
the r to four pions. The model is based on a fit to e + e~ data from Novosibirsk. 

32 It should be noted that there were a number of mistakes in this paper, which were corrected in Ref. [118]. 
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FivePionCurrent We use the model of Ref. [126], which includes pu and pa intermediate 
states, via the eti meson to model the five-pion decay modes of the r. 

9.4 Strong and electromagnetic hadron decays 

The vast majority of the strong and electromagnetic decays in Herwig++ are simulated using a few 
simple models based on the spin structure of the decay. These simple models are supplemented 
with a small number of specialized models, usually from experimental fits, for specific decay 
modes. In this section we describe the different models we use for these decays for the scalar, 
vector and tensor mesons. All of these are implemented in Decayer classes that inherit from the 
Decayl integrator class of Herwig++. 

For a number of the decays of bottomonium and charmonium resonances we use inclusive 
electromagnetic and strong decays to qq, gg, ggg and ggj, which are described in a separate 
section. 

A number of decays are still performed using a phase-space distribution generated using the 
Hw64Decayer, which implements the same models as were available in the FORTRAN HERWIG 
program. In addition we use the MAM BO algorithm, [127], implemented in the MamboDecayer 
class, to generate the momenta of the decay products according to a phase-space distribution for 
a number of high-multiplicity modes. 

9.4.1 Scalar mesons 

While the majority of the scalar meson decays are performed using general Decayers based on 
the spin structures there are a number of models implemented for the rare radiative decays of 
the light pseudoscalar mesons, three-body decays of the rj and rj', and the decay ir° — > e + e~e + e~ . 

EtaPiGammaGammaDecayer We use the Vector-Meson Dominance (VMD)-based model 
of Ref. [128] for the decays rj, rj' — > 7r°77. In practice we use a running width for the p to include 
the rj' decay as well as the rj decay and take the parameters from Ref. [128]. 

EtaPiPiGammaDecayer We use either a VMD type model or a model using either the the- 
oretical or experimental form of the Omnes function 33 taken from Refs. [128, 129] for the decay 
of the rj or rj' to 7r + 7r~7. 

EtaPiPiPiDecayer The decay of a pseudoscalar meson, for example the rj or rj', to two charged 
and one neutral or three neutral pions, or of the rj' to two pions and the rj, is performed using a 
parameterization of the matrix element squared taken from Ref. [130]. The experimental results 
of Refs. [131] and [132] are used for the rj — > 7r + 7r~7r° and rj — ► 7r 7r 7r° decays respectively. The 
theoretical values from Ref. [130] are used for the other decays. 

PScalar4FermionsDecayer As the ir° is so copiously produced it is one of the small number 
of particles for which we include branching ratios below the level of 10~ 4 . The matrix element 
for the sub-leading decay 7r° — > e + e~e + e~ is taken to be the combination of the standard matrix 
element for 7r° — > 77 and the branching of the photons into e + e~ . 

33 0ur default choice is to use the experimental form of the Omnes function. 
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PScalarP Scalar VectorDecayer This matrix element is used to simulate the decay of the 2S 
pseudoscalar mesons to a vector meson and a IS pseudoscalar meson. It is also used for the 
decay of some scalar mesons to vector mesons and another scalar meson, which has the same 
spin structure. The matrix element has the form 

M = ge^ Po + Pl )^ (9.11) 

where e 2 is the polarization vector of the vector meson, po is the momentum of the decaying 
particle, p\ is the momentum of the outgoing pseudoscalar meson and g is the coupling for the 
decay. 

P Scalar VectorFermionsDecayer There are a number of decays of a pseudoscalar meson to 
either a vector meson or the photon and a lepton-antilepton pair. The classic example is the 
Dalitz decay of the neutral pion, 7r° — > / ye + e~ . We take the propagator of the off-shell photon to 
be — ij-, where m f j is the mass of the fermion-antifermion pair. The option of including a vector 

meson dominance form factor is included. 

P Scalar Vector VectorDecayer In practice the vast majority of the decays of pseudoscalar 
mesons to two spin-1 particles are of the form P — > 77 for which, because the photon is stable, 
it is not as important to have a good description of the matrix element. There are however some 
decays, e.g. rj' — > CJ7, for which this matrix element is needed. 
The matrix element is taken to be 

M = ge^ Plll e lu p 2a e 2f3 , (9.12) 

where pi t2 and ei )2 are the momenta and polarization vectors of the outgoing vector particles and 
g is the coupling for the decay. 

ScalarMesonTensorScalarDecayer There are a limited number of decays of a 
(pseudo) scalar meson to a tensor meson and another (pseudo) scalar meson. The matrix ele- 
ment takes the form 

M=ge a ^p 0a p 2 ^ (9.13) 

where e a/3 is the polarization tensor of the outgoing tensor meson, p is the momentum of the 
decaying particle, p 2 is the momentum of the outgoing (pseudo) scalar meson and g is the coupling 
for the decay. 

ScalarScalarScalarDecayer The decay of a scalar meson to two scalar mesons has no spin 
structure and we assume that the matrix element is simply constant, i.e. 

M = g. (9.14) 

We still include a matrix element for this decay in order to simulate both the off-shell effects in 
the decay and to give the correct partial width to be used in the running width calculation for 
the incoming particle. 
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Scalar VectorVectorDecayer A number of the scalar mesons decay to two vector mesons. 
The matrix element is taken to have the form 

M = g [pi ■ p 2 e x -e 2 -pi- e 2 p 2 ■ e x ] , (9.15) 

where e^ 2 are the polarization vectors of the outgoing vector particles and pi >2 are their momenta. 

9.4.2 Vector mesons 

With the exception of the three-pion decay modes of the u, <fr and a\ mesons, and the two-pion 
decays of onium resonances, we use general Decayers based on the spin structure for all the strong 
and electromagnetic vector and pseudovector meson decays. 

alSimpleDecayer This class implements the model of Kiihn and Santamaria [119] for the 
decay of the a\ meson to three pions and only includes the lightest two p meson multiplets in 
the modelling of the decay. 

alThreePionCLEODecayer This class implements the model of CLEO [123] for a\ decay to 
three pions, which is a fit to CLEO data on r~ — > 7r 7r 7r~iv The model includes the coupling 
of the a\ to the p, p(1450), / (1370), / 2 (1270) and a mesons. 

alThreePionDecayer This class implements a model of a± decay to three pions based on the 
modelling of the a\ used in the 47r currents for tau decays presented in Ref. [125] and includes 
the p and a resonances. 

OniumToOniumPiPiDecayer The decay of onium resonances to lighter states and a pion 
pair, O' — > Otttt, uses the matrix element [133] 

M = e' ■ e [A (q 2 - 2m\) + BE 1 E 2 ] + C ((e' • 9l )(e • q 2 ) + (e' • q 2 )(e ■ Ql )) , (9.16) 

where e' is the polarization vector of the decaying onium resonance, e is the polarization vector 
of the outgoing onium resonance, A, B and C are complex couplings, is the pion mass, E± j2 
are the pion energies, qi t2 are the pion momenta and q is the momentum of the ttti system. 

The results of BES [134] are used for $ -> J/i/j and CLEO [135] for T(3S) and T(2S) 
decays. The remaining parameters are chosen to approximately reproduce the distributions from 
BaBar [136] and CLEO [137] for T(4S) and ^(3770) decays respectively. 

PVectorMesonVectorPScalarDecayer The matrix element for the decay of a pseudovector 
meson to a spin-1 particle, either a vector meson or a photon, and a pseudoscalar meson is taken 
to be 

M = ge„ [p v ■ p ey - Pv e v • Po) , (9.17) 

where ey is the polarization vector of the outgoing vector meson, py is the momentum of the 
outgoing vector meson, e is the polarization vector of the decaying pseudovector and po is the 
momentum of the decaying particle. 
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VectorMeson2FermionDecayer Most of the decays of the vector mesons to a fermion- 
antifermion pair are the decays of the light vector mesons to electron and muon pairs, and 
of the bottomonium and charmonium resonances to all the charged leptons. In addition we use 
this matrix element for some baryonic charmonium decays. 
The matrix element is taken to have the form 



M=ge tt u(p f )'fv(p f ) 1 (9.18) 

where g is the coupling for the decay, p/ is the four-momentum of the outgoing fermion, pj is 
the four-momentum of the outgoing antifermion and e is the polarization vector of the decaying 
particle. 



VectorMeson2MesonDecayer The matrix element for the decay of a vector meson to two 
scalar (or pseudoscalar) mesons is given by 

M = gyppe ■ (p! - pa), (9-19) 

where gypp is a dimensionless coupling, e is the polarization vector of the decaying particle and 
P12 are the momenta of the outgoing scalars. 



VectorMeson3PionDecayer Both the lowest-lying isospin-zero vector mesons, u and 0, have 
large branching ratios for the decay into three pions. For these mesons the decay is assumed 
to be dominated by the production of the lowest lying p multiplet. Our default model for the 
matrix element for this decay is 



M 



ge^ 131 'e^p la p 2 pP3„ 



d + ^h [BW^i 2 ) + BW<(si3) + BW,( S23 )] 



(9.20) 



where pi,2,3 are the momenta of the outgoing particles, s^- = (p» +pj) 2 , g is the overall coupling 
for the decay, d is a complex coupling for the direct interaction, is the coupling of the zth p 
multiplet and BWj(s) is a Breit-Wigner distribution with a p-wave running width. This is an 
extension of the model used by KLOE [138] to include higher p multiplets. 



VectorMesonP Scalar FermionsDecayer The decay of a vector meson to a pseudoscalar me- 
son and a fermion- antifermion pair is simulated using a matrix element based on that for the 
V — > VP vertex combined with the branching of the vector, which is in reality always a photon, 
into a fermion-antifermion pair. 



VectorMesonP VectorPScalarDecayer There are a number of decays of both the charmo- 
nium resonances and light vector mesons from the higher multiplets to pseudovector mesons. 
The matrix element for the decay is 

M = g [pa ■ Po^a -e-p A -ee A - p ] , (9.21) 



where e A is the polarization vector of the outgoing pseudovector meson, p A is its momentum, e 
is the polarization vector of the decaying particle and po is its momentum. 
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VectorMesonVectorPScalarDecayer The decay of a vector meson to another spin-1 particle 
and a pseudoscalar meson is common in both the radiative decay of the IS vector mesons and 
the decay of higher vector multiplets to the IS vector mesons. The matrix element for the decay 
is 

M = ge^ a(S e 0f ,p 0u p la e 1(3 , (9.22) 

where po is the momentum of the decaying particle, p\ is the momentum of the outgoing vector 
particle, e is the polarization vector of the incoming meson and €\ is the polarization vector of 
the outgoing vector particle. 



VectorMesonVectorScalarDecayer We include a number of decays of the vector mesons to 
a scalar meson and either the photon or another vector meson. In practice the vast majority 
of these decays are radiative decays involving scalar mesons. The remaining decays use the a 
meson as a model for four-pion decays of the excited p multiplets. 
The matrix element for the decay is 

M = ge„ [p v ■ p e^ - p^e v ■ p ) , (9.23) 

where g is the coupling for the decay, e is the polarization vector of the decaying vector meson, 
ey is the polarization vector of the outgoing vector meson, p is the momentum of the decaying 
particle and py is the momentum of the outgoing vector meson. 



VectorMesonVectorVectorDecayer There are a small number of decays of excited p mul- 
tiplets to p mesons included in the simulation. We model these decays using the matrix element 



M 



g 



a 



(9.24) 



where g is the coupling for the decay, ei^ are the polarization vectors of the outgoing mesons, 
Pi t 2 are the momenta of the outgoing mesons, e is the momentum of the decaying particle and po 
is its momentum. 



9.4.3 Tensor mesons 

Only a relatively small number of tensor meson states are included in the simulation, compared 
to the vector and scalar mesons. All their decays are simulated using a small number of matrix 
elements based on the spin structure of the decays. Many of the multi-body decays of the tensor 
mesons are simulated using these two-body matrix elements with off-shell vector and scalar 
mesons. 



TensorMeson2PScalarDecayer The matrix element for the decay of a tensor meson to two 
pseudoscalar (or scalar) mesons is 

M = ge^p!^, (9.25) 

where g is the coupling for the decay, pi^ are the momenta of the decay products and e^ u is the 
polarization tensor for the decaying meson. 
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TensorMesonVectorPScalarDecayer There are a number of decays of tensor mesons to a 
spin-1 particle, either a vector meson or the photon, and a pseudoscalar meson, examples include 
a 2 — > pir and a 2 — > ir'-f. The matrix element is taken to be 

M = e^pp^afrp^pl, (9.26) 

where g is the coupling for the decay, pp is the momentum of the pseudoscalar meson, py is the 
momentum of the vector, ey is the polarization vector of the outgoing vector meson and e^ u is 
the polarization tensor for the decaying meson. 



TensorMesonVectorVectorDecayer We have based our matrix element for the decay of a 
tensor meson to two vector mesons on the perturbative graviton decay matrix element [139] in 
such a way that it vanishes if the polarizations of the outgoing vectors are replaced with their 
momenta. The matrix element is 



M 



+ (A* ~ *)} 



(9.27) 



a P 
e 2P2 



P a 
^2 



where g is the coupling for the decay, ex s 2 are the polarization vectors for the outgoing vector 
mesons and e^ u is the polarization tensor for the decaying meson. In practice this matrix element 
is mainly used with off-shell vector mesons to model three- and four-body decays of the tensor 
mesons. 



9.4.4 Baryon Decays 

The strong and electromagnetic decays of the baryons are modelled in Herwig++ using models 
based on either heavy quark effective theory, for the baryons containing a bottom or charm quark, 
or SU(3) symmetry for the light baryons. 

All the strong decays, and many of the weak hadronic decays, involve the decay of a spin-| 
or | baryon to either a pseudoscalar meson or a vector particle and another spin- ^ or | baryon. 
Lorentz invariance restricts the possible form of the matrix elements. We use the following matrix 
elements which are implemented in the BaryonlMesonDecayerBase class from which the Decayers 
inherit. We use the following matrix elements 



M 


= u{p l ){A + B l5 )u(p ) 


1 

2 ~ 


1 

♦ 5 + o 


(9.28a) 


M 


= u( Pl )e*P [ 7/9 (Ai + B xlb ) + p p(A 2 + B 2l5 )} u(p ) 


1 

2 ~ 


1 


(9.28b) 


M 


= ^(p^poa [A + .675] u(p ) 


1 

2 ~ 


♦ 5 + o 


(9.28c) 


M 


= u a ( Pi yP [g a(3 (A 1 + B ll5 ) 


1 

2 ~ 


4 + 1 


(9.28d) 



+P0a(A 2 + £275) + POaPop(A 3 + £375)] U(p ) 
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for spin- 2 decays and 

M = U(pi)pia [A + ,675] U a ( Po ) 

M = uipje* 13 \g aP {A 1 + Si 7B ) 

+Pi a (A 2 + .6275) + Pi a Pw(Az + 5375)] u a (p ) 
M = u a ( Pl ) [(A! + B ll5 )g a p + p 0a Pip{A 2 + B 2l5 )} u^po) 

for spin-| decays. In general u(po) is the spinor of a decaying spin-| baryon, u 13 ^) is the spinor 
of a decaying spin-| baryon, u{pi) is the spinor for an outgoing spin-| baryon and u^(px) is 
the spinor for an outgoing spin-| baryon. The momentum of the decaying baryon is po, of the 
outgoing baryon is p\ and of the outgoing meson is p 2 - All the matrix elements are parameterized 
in terms of a number of coefficients A and B which can in principle depend on the momentum 
transfered in the decay. 

RadiativeHeavyBaryonDecayer There are a number of radiative decays of heavy baryons 
included in the simulation. Apart from some transitions of charm baryons, e.g. "El c — > H c 7, these 
transitions have not been observed and are included based on model calculations based on heavy 
quark effective theory [115]. 

RadiativeHyperonDecayer The radiative decays of hyperons are simulated using the model 
of Ref. [140]. 

StrongHeavyBaryonDecayer The StrongHeavyBaryonDecayer class implements the strong 
decays of bottom and charm baryons using the results of Ref. [115]. 

SU3BaryonDecupletOctetPhotonDecayer The decay of a decuplet baryon to an octet 
baryon and a photon is assumed to occur via the SU(3) conserving term in the chiral Lagrangian. 

SU3BaryonDecupletOctetScalarDecayer This Decayer is based on SU(3) symmetry for 
the decay of a decuplet baryon to an octet baryon and a scalar meson. 

SU3BaryonOctetDecupletScalarDecayer The SU3BaryonOctetDecupletScalarDecayer per- 
forms the decay of excited octet baryons to decuplet baryons and a scalar meson using a La- 
grangian based on SU(3) symmetry. 

SU3BaryonOctetOctetPhotonDecayer The SU3BaryonOctetOctetPhoton Decayer models 
the radiative decay of excited octet baryons using a Lagrangian based on SU(3) symmetry. 

SU3BaryonOctetOctetScalarDecayer The SU3BaryonOctetOctetScalarDecayer simulates 
the strong decay of excited octet baryons using a Lagrangian based on SU(3) symmetry. 



3 1 . . 

- - - + (9.28e) 

+ ( 9 -28g) 
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SU3BaryonSingletOctetPhotonDecayer The SU3BaryonSingletOctetPhotonDecayer models 
the radiative decay of singlet baryons using a Lagrangian based on SU(3) symmetry. 

SU3BaryonSingletOctetScalarDecayer The SU3BaryonSinglet0ctetScalarDecayer simulates 
the strong decay of singlet baryons using a Lagrangian based on SU(3) symmetry. 

9.4.5 Inclusive strong and electromagnetic decays 

For a number of bottomonium and charmonium resonances we make use of partonic decays of 
the mesons to model the unobserved inclusive modes needed to saturate the branching ratios. 
These decays are modelled using the QuarkoniumDecayer class, which implements the decay of 
the onium resonances to qq and gg according to a phase-space distribution, and the decay to ggg 
and ggj according to the Ore-Powell matrix element [141]. The QuarkoniumDecayer class inherits 
from the PartonicDecayerBase, which uses the cluster model to hadronize the resulting partonic 
final state with a veto to ensure that there is no double counting with the exclusive modes. 

9.5 Weak hadronic decays 

There are five classes of weak mesonic decays currently included in the simulation: 

1. weak exclusive semi-leptonic decays of bottom and charm mesons; 

2. weak exclusive hadronic decays of bottom and charm mesons; 

3. weak inclusive decays; 

4. weak leptonic decay of pseudoscalar mesons; 

5. weak inclusive b — > sj mediated decays. 

We adopt a number of different approaches for these decays as described below. 

9.5.1 Exclusive semi-leptonic decays 

The matrix element for exclusive semi-leptonic decays of heavy mesons, X Y£u, can be written 

as 



where p^ is the momentum of the outgoing charged lepton, p v is the momentum of the neutrino 
and Gp is the Fermi constant. The hadronic current {X\{V — A)^\Y) can be written as a general 
Lorentz structure, for a particular type of decay, with a number of unknown form factors. 

We have implemented a number of form-factor models based on experimental fits, QCD sum 
rule calculations and quark models. The form factors for the weak decay of pseudoscalar mesons 
are implemented using the general Lorentz-invariant form. In each case the momentum of the 
decaying particle, X, is px while the momentum of the decay product, Y, is py. In general the 
form factors are functions of the momentum transfer q 2 where q = Px ~ Vy ■ The masses of the 
decaying particle and hadronic decay product are mx and my respectively. 




(9.29) 



104 



The scalar-scalar transition matrix element is defined by 



2 2 

m v — mi 



(Y(p Y )\(V - A),\X{p x )) = ( P x+Py),U(Q 2 ) + | X g2 Y jq, [fo(q 2 ) ~ U{q 2 )] , (9.30) 

where f+(q 2 ) and fo(q 2 ) are the form factors for the transition. In general the terms proportional 
to g M give rise to contributions proportional to the lepton mass for semi-leptonic decays and 
therefore only contribute to the production of tau leptons. 
The scalar-vector transition matrix element is defined to be 

(Y(p Y )\(V-A),\X(p x )) = (9.31) 

A 2 {q 2 ) 



ie^mx + m Y )Ax(q ) + i(p x + Py)^* ■ q 



m x + my 



+iq fl e ■ q^- (A 3 (q ) - A Q (q )) + e^ pu e p' x p Y - 



q z m x + my 

where the form factor A^{q 2 ) can be defined in terms of A\ and A-i using 

2 m x + my 2 m x - m y 2 n / n qo n 

Mq ) = ) - ~^-Mq ) 0.32) 

and A (0) = A 3 (0). The independent form factors are A (q 2 ), Ai(q 2 ), A 2 {q 2 ) and V(q 2 ). 
The transition matrix element for the scalar-tensor transition is 



(Y(p Y )\(V-A) ll \X(p x ))= (9.33) 

* 



th(q 2 )e, vXp e*» a p Ya (p x + P y) A (p x - Vy) P ~ Hq^Py 
-h(q 2 y a pPxPx(Px + Py)„ ~ b.{q 2 )e* a(3 p x p x {p x - p Y ) 



where h(q 2 ), k(q 2 ), b_(q 2 ) and b + (q 2 ) are the Lorentz invariant form factors. 

The combination of the form factors and the leptonic current is handled by the 
SemiLeptonicScalarDecayer class, which combines the form factor and the current to calculate 
the matrix element and uses the methods available in the Decay Integrator class, from which it 
inherits, to generate the momenta of the decay products. 

In addition to the form factors for the standard weak current we include the form factors 
needed for weak radiative decays where available, although these are not currently used in the 
simulation. 

The various form factors that are implemented in Herwig++ are described below. They all 
inherit from the ScalarFormFactor class and implement the relevant virtual member functions for 
the calculation of the form factors. 

BallZwickyScalarFormFactor This is the implementation of the QCD sum rule calculation 
of the form factors of Ref. [142] for the decay of a S-meson to a light pseudoscalar meson. 

BallZwicky VectorFormFactor This is the implementation of the QCD sum rule calculation 
of the form factors of Ref. [143] for the decay of a 5-meson to a light vector meson. 
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HQETFormFactor This implements the form factors for the transitions between mesons con- 
taining bottom and charm quarks in the heavy quark limit. The parameterization of Ref. [144] 
for the finite-mass corrections is used together with the experimental results of Refs. [145, 146]. 

ISGWFormFactor The ISGW form factor model [147] is one of the original quark models 
for the form factors and is included in the simulation mainly for comparison with the later, 
ISGW2 [148], update of this model. This set of form factors has the advantage that it includes 
form factors for most of the transitions required in the simulation. The form factors are taken 
from Ref. [147] together with the form factors that are suppressed by the lepton mass from 
Refs. [149,150]. 

ISGW2FormFactor The ISGW2 form factors [148] are an update of the original ISGW form 
factors [147]. As with the original model they are based on a quark model and supply most of 
the form factors we need for the simulation. 

KiselevBcFormFactor This is the implementation of the form factors of Ref. [151] for the 
weak decays of B c mesons. This model is used as the default model for weak B c decays as the 
branching ratios for the B c meson used in the simulation are calculated using the same model. 

MelikhovFormFactor This is the implementation of the relativistic quark model form factors 
of Ref. [152] for B -> n,p. 

MelikhovStechFormFactor This is the implementation of the model of Ref. [153], which is 
an update of the model of Ref. [152] including additional form factors. 

WSBFormFactor This is the implementation of the form factor model of Ref. [154] for the 
semi-leptonic form factors. It includes form factors for a number of D, B and D s decays. In 
practice the parameters of the model were taken from Ref. [155], which includes a number of 
transitions that were not considered in the original paper. 

This form factor model is included both to give an alternative for many modes to the ISGW 
models and for use in the factorization approximation for hadronic charm meson decays. 

We also include exclusive semi-leptonic decays of heavy baryons in the same way. The tran- 
sition matrix elements are given below for the decay X{p x ) Y(Py) with = (p x — Py) 1 mj as 
for the mesonic case. The transition matrix for the | — ► | transition is defined as 



(Y(p Y )\(V-A),\X(p x )) 



i 



(m + mi) 



1 



(t {FX + F 3 A l5 ) u( Px ) 



(9.34) 



(m + mi) 
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The transition matrix element for the ~ — > | transition is 



(Y(p Y )\{V-A)»\X(p x )) = (9.35) 
u a ( P y) 



g afl (G\ + Gf l5 ) + {mo l mi) Po a l, (G v 2 + Gfr*) 



where again the dependence of the form factors GV 2 3 4 on the momentum transfer q 2 has been 
suppressed. These definitions differ from those in the literature because we have divided some 
terms by the sum of the baryon masses to ensure that the form-factors are all dimensionless. 

We have implemented a number of different models for the baryon form factors, mainly based 
on quark model calculations. All the form factor classes inherit from the BaryonFormFactor class 
and implement the calculation of the form factors in the specific model. The 
SemiLeptonicBaryonDecayer class handles the combination of the form factor and the leptonic 
current to calculate the partial width and decay kinematics. 

The models we have implemented are: 

BaryonSimpleFormFactor This is a simple form factor model for for the semi-leptonic decay 
of the light baryons. The form factors are assumed to be constant and are taken from the quark 
model results of [156]. 

BaryonThreeQuarkModelFormFactor This model is the implementation of the relativistic 
three-quark model calculation of [109] for the form factors of baryons containing a heavy quark. 

As the only formulae in the paper, in a form which can be implemented in the simulation, 
are for the heavy-to- heavy i.e. bottom to charm decays these are the only modes included, 
although the paper also includes charm decays and bottom decays to light quarks. The form 
factors are calculated by numerical computing the integrals from [109] to obtain the coefficients 
for an expansion of the form factors in uj. 

ChengHeavyBaryonFormFactor This is a quark model calculation [114,157] of form factors 
for bottom and charm baryons. It is used for some bottom and charm baryon semi-leptonic 
decays. However it is mainly intended to implement the factorization approximation results 
of [114] for non- leptonic decays. 

LambdabExcitedLambdacSumRuleFormFactor This is the QCD sum rule based calcu- 
lation of [113] for the form factors for the decay of the A° to excited A+ states. This is used 
for the semi-leptonic decay of the the A° to excited A+ states to model the part of the total 
semi-leptonic branching ratio of the A° not accounted for by the production of the A+. 

LightBaryonQuarkModelFormFactor This is a relativistic quark model calculation [158] 
of the form factors for the decay of baryons containing the light quarks. 

SingletonFormFactor This model is a quark model calculation [159] of the form factors of 
spin-1 baryons containing a bottom or charm quark. 
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9.5.2 Exclusive hadronic decays 

We include two types of simulation of exclusive weak meson decays. The first is based on the 
naive factorization approximation. If we consider, for example, the decay of a charm meson then 
the effective Lagrangian for the interaction can be written as [155] 

Gp 

£ c fr = -^VudVsc K {uj^P L d) (s'J^Plc) + a 2 (f7 M P £ d) {uj^Plc)) , (9.36) 

where Gp is the Fermi constant, V u d and V sc are the relevant CKM matrix elements and 
are scale-dependent coefficients. The remainder of the expression involves the currents for the 
quark fields. When we consider the transition between mesonic states the matrix element can be 
written in terms of the form factors, for the c — > s or c — > u transitions, and weak currents for 
(u'j^PLd) or (s-i^P L d). 

This allows us to simulate weak hadronic decays using the form factors we have already im- 
plemented for semi-leptonic meson decays together with the weak currents from tau decays. The 
combination of the form factor classes, which inherit from ScalarFormFactor, and weak currents, 
which inherit from WeakDecayCurrent, is handled by the ScalarMesonFactorizedDecayer class for 
the simulation of hadronic weak meson decays in the factorization approximation. Similarly the 
combination of weak form factors inheriting from the BaryonFormFactor class and weak currents is 
handled by the BaryonFactorizedDecayer class for the simulation of hadronic weak baryon baryon 
decays in the factorization approximation. 

In addition to the weak exclusive decays based on the factorization approximation we include 
a small number of classes for the simulation of D — > Ktttt Dalitz decays based on various 
experimental fits. Currently there are three such models implemented. 

DtoKPiPiCLEO This class implements the CLEO fits of Refs. [160] and [161] for the decays 
D° — > K°TT + ir~ and D° — > K~7t + tt°. This is our default simulation of these decays. 

DtoKPiPiE691 The DtoKPiPiE691 class implements the model of E691 [162] for the decays 
D° — > K°7r + 7T~, D° — > K~7t + 7t° and D + — > Ji~7r + 7r~. This is our default simulation for the 
D + — > K'-k + -k~ decay. 

DtoKPiPiMarklll This class implements the model of the Mark-Ill collaboration for the 
decays D° -> K°n + n~, D° -> K~tt+tt°, D+ -> R-iT+it- and D+ -> K°n + Tr°. This is our default 
model for the decay mode D + — > K°7t + ii . 

While some of the exclusive weak hadronic decays are simulated using the factorization ap- 
proximation we also use a number of other models which include non-factorizable contributions. 
These all inherit from Baryon 1 Meson DecayerBase which performs the calculation of the matrix 
elements. 

KornerKramerCharmDecayer This is the implementation of the results of the spectator 
quark model of [108] for the non-leptonic weak decay of charm baryon, i.e. A)C + , H°, and 

This model provides branching ratios and decay matrix elements for a large number of modes 
and is used as the default simulation for many of the hadronic decay modes of the weakly 
decaying charm baryons. In addition as for many of these baryons all the decay modes have not 
been observed in some cases the branching ratio calculations are used to add these modes. 
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NonLeptonicHyperonDecayer We use the results of [163] for the matrix elements for the 
weak, non-leptonic, decay of a number of hyperons, i.e. S ± '°, S 0, ~. 

The matrix element for the decay is given in terms of the invariant amplitudes 

C = u Bj {A + B l5 }u Bi (9.37) 

where Bj is the outgoing baryon and Bi is the incoming baryon. 
The default amplitudes are taken from the fit in [163]. 

NonLeptonicOmegaDecayer We use the model of [164] for the non-leptonic weak decays of 
the Or to a baryon from the lightest SU(3) octet and a pseudoscalar meson. Due to problems 
with the size of the d-wave term in this model, and recent measurements giving the opposite sign 
for the a parameter, we have set this term to zero in the simulation. 



OmegaXiStarPionDecayer We use the model of [165] for the weak decay of the Q~ to the 
H* and a pion. This decay has a very low branching ratio and the model is mainly included to 
test the code involving the decay of a spin-| particle to another spin-| particle. 

9.5.3 Weak inclusive decays 

In addition to the exclusive weak decays of the hadrons to specific final states we include a number 
of models of the decay of hadrons containing a heavy, i.e. bottom or charm, quark based on the 
partonic decay of the heavy quark. The Herwig++ cluster hadronization model is then applied 
to the resulting partonic final state to produce hadrons. This approach is primarily used for the 
bottom hadrons where there are insufficient exclusive modes to saturate the branching ratios. 
All of the classes implementing partonic decay models inherit from the PartonicDecayerBase to 
use the cluster hadronization model to hadronize the partonic final state. 

The HeavyDecayer class implements the weak decays of hadrons using either the weak V — A 
matrix element or flat phase space. The WeakPartonicDecayer includes additional features to 
simulate decays intended to increase the rate of baryon production and gluonic penguin decays. 

In addition the BtoSGammaDecayer for weak penguin-mediated decays, described in Sect. 9.5.5, 
and the QuarkoniumDecayer class for the decay of bottomonium and charmonium resonances, de- 
scribed in Sect. 9.4.5, also perform partonic decays and inherit from the PartonicDecayerBase 
class. 



9.5.4 Leptonic decays 

There are a small number of decays of pseudoscalar mesons to a charged lepton and a neutrino, 
e.g. 7r — > Iv and D s — > iv. For most of these decays the inclusion of the matrix element is 
superfluous as the decay products are stable. However the B and D s mesons can decay in this 
way to the r and therefore we include the PScalarLeptonNeutrinoDecayer class to simulate these 
decays using the matrix element 

M = -j=jpG F V CK Mrriiu{pi){l - ls)v{p u ), (9.38) 

where fp is the pseudoscalar decay constant, Gp is the Fermi constant, Vckm is the relevant 
CKM matrix element, is the mass of the lepton, is the momentum of the charged lepton 
and p u is the momentum of the neutrino. 
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9.5.5 b -> s 7 

There is a range of decays, both inclusive and exclusive, mediated by the b — > 57 transition. 
We currently only include modelling of the inclusive decay. These decays are simulated by using 
a partonic decay of the B meson to a photon and a hadronic system, composed of a quark 
and antiquark, which recoils against the photon. The mass spectrum of the hadronic system is 
calculated using a theoretical model. 

The calculation of the mass spectrum is handled by a class inheriting from the 
BtoSGammaHadronicMass class. Different models of the mass spectrum can then be implemented 
by inheriting from this class. Currently we have only implemented two such models. The first, 
BtoSGammaFlatEnergy, is solely designed for testing and generates a mass spectrum such that 
the photon energy distribution is flat. The second model, BtoSGammaKagan, which is the de- 
fault, implements the theoretical calculation of Ref. [166]. The BtoSGammaDecayer then uses 
the calculation of the hadronic mass spectrum to simulate the partonic decay as a model of the 
inclusive mode. As with the Decayers described in Sect. 9.5.3 the BtoSGammaDecayer inherits 
from the PartonicDecayerBase class to use the cluster model to perform the hadronization of the 
partonic final state. 

9.6 Code structure 

The HwDecayHandler, which inherits from the DecayHandler class of ThePEG, is responsible for 
handling all particle decays in Herwig++. It uses the DecaySelector from the ParticleData object 
of the decaying particle to select a DecayMode object corresponding to a specific decay according 
to the probabilities given by the branching ratios for the different modes. The DecayMode object 
then specifies a Decayer object that is responsible for generating the kinematics of the decay 
products for a specific decay. 

All of the Decayer classes in Herwig++ inherit from the HwDecayerBase class, which in turn 
inherits from the Decayer class of ThePEG. In turn, with the exception of the Hw64Decayer and 
MamboDecayer classes, which implement general phase-space distributions for the decay products, 
all the Decayer classes in Herwig++ inherit from either the Decay Integrator or PartonicDecayBase 
classes. 

The Decayl integrator class provides a sophisticated multi-channel phase space integrator to 
perform the integration over the phase space for the decays. This means that the calculation 
of the matrix element and specification of the phase-space channels are all that is required to 
implement a new decay model. The majority of the matrix elements are calculated as helicity 
amplitudes, which allows the spin-propagation algorithm of Refs. [27-30] to be implemented. 
The structure of the Herwig++ Decayer classes and HwDecayHandler is designed so that these 
correlations are automatically included provided the helicity amplitudes for the matrix elements 
are supplied. 

The PartonicDecayBase class provides a structure so that the decay products of a partonic 
hadron decay can be hadronized using the cluster model, while at the same time ensuring that 
there is no overlap with the particle's exclusive decay modes. All classes implementing partonic 
decays in Herwig++ inherit from the PartonicDecayBase class. 

Certain Decayer classes also make use of helper classes to implement the decays. The main 
examples are: 
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• the WeakDecayCurrent class, which provides a base class for the implementation of weak 
hadronic currents, is used by the TauDecayer, SemiLeptonicScalarDecayer, 
SemiLeptonicBaryonDecayer, ScalarMesonFactorizedDecayer and BaryonFactorizedDecayer 

classes, which implement tau decays, semi-leptonic meson and baryon decays and hadronic 
weak meson and baryon decays using the naive factorization approximation, respectively; 

• the Sea larForm Factor class, which provides a base class for the implementation of the scalar 
form factors and is used by the SemiLeptonicScalarDecayer and ScalarMesonFactorizedDecayer 
classes, which implement semi-leptonic meson decays and hadronic weak meson decays 
using the naive factorization approximation, respectively; 

• the BaryonFormFactor class, which provides a base class for the implementation of the 
baryon form factors and is used by the SemiLeptonicBaryonDecayer and 
BaryonFactorizedDecayer classes, which implement semi-leptonic baryon decays and hadronic 
weak baryon decays using the naive factorization approximation, respectively; 

• the BtoSGammaHadronicMass, which provides a model of the hadronic mass spectrum in 
inclusive b —>■ 37 decays performed by the BtoSGammaDecayer class. 

The vast majority of the decay models have a large number of parameters, all of which 
are accessible via the Interfaces of the classes. A more detailed description of both the physics 
models used in the code and their parameters can be found in the Doxygen documentation and 
Refs. [25,92]. 

There are a number of classes that are designed to include the off-shell weight given in 
Eq. (9.1) in the generation of the particle decays. The GenericWidthGenerator is designed to 
use the information on the partial widths for the different decay modes supplied by the Decayer 
classes, which inherit from Decay Integrator, to calculate the running width for a given particle. 
The GenericMassGenerator class then uses the running width to allow the weight given in Eq. (9.1) 
to be included when generating the particle decays. The inheriting ScalarMassGenerator class 
implements the Flatte lineshape [116] for the ao(980) and /o(980) mesons. 

For decays where the decay products can be off-shell, and three-body decays, integrals over 
either the masses of the decay products or the three-body phase space must be performed in 
order to calculate the running partial widths. In order to facilitate the calculation of the partial 
widths a number of classes inheriting from the WidthCalculatorBase class are implemented to 
calculate the partial widths for various decays: 

• the TwoBodyAIIOnCalculator returns the partial width for a two-body decay where both the 
decay products are on mass-shell; 

• the OneOffShellCalculator returns the partial width for a decay where one of the outgoing 
particles is off mass-shell; 

• the TwoOffShellCalculator returns the partial width for a decay where two of the outgoing 
particles are off mass-shell; 

• the ThreeBodyAIIOnCalculator returns the partial width for a three-body decay where all 
the decay products are on mass-shell by performing the two non-trivial integrals over the 
phase-space variables; 
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• the ThreeBodyAIIOnllntegralCalculator returns the partial width for a three-body decay 
where all the decay products are on mass-shell by performing one integral over the phase- 
space variables, this requires that the second integral has already been performed analyti- 
cally. 
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10 Summary 



In this manual we have described the physics and structure of Herwig++ version 2.3. More 
detailed technical documentation can be obtained from the web site 

http : // pro j ects . hepf orge . org/herwig 

as well as a growing number of user guides, example applications, frequently-asked-questions and 
other useful information. Most of this is obtained by following the "wiki" link at the top of the 
page. To be able to contribute to the wiki and submit trac tickets, please email the authors, at 

herwig@proj ects .hepf orge . org 

To improve the current version of Herwig++ and plan development of future versions, we depend 
on feedback from users. If you use Herwig++ please register at the address above and post your 
experience (positive or negative) and code examples you feel other users would benefit from, and 
open a trac ticket for any bugs or unexpected features you find, as well as any new features or 
improvements you would like to see. Of course, for any bug report, the more clearly you can 
illustrate the problem, and the fact that it is a problem with Herwig++ and not an external 
package it is connected to, the more quickly we are likely to be able to solve it. 

Herwig++ has been extended enormously since the last version for which a published manual 
exists, 1.0. It now provides complete simulation of hadron-hadron collisions with a new coherent 
branching parton shower algorithm, including quark mass effects, a sophisticated treatment of 
BSM interactions and new particle production and decay, an eikonal model for multiple partonic 
scattering, greatly improved secondary decays of hadrons and tau leptons and a set of input 
parameters that describe e + e~ annihilation data rather well. 

New features planned for the near future include: an improved treatment of baryon decays; 
spin correlations within the parton shower; 'multiscale' showering of unstable particles; simulation 
of DIS processes; B mixing; and an improved treatment of gluon splitting to heavy quarks. Of 
course we are all users of Herwig++ as well as developers and are working on a large number 
of other new features related to phenemenological studies we are making. The list will continue 
to grow, according to the physics interest and needs of ourselves and others using it for physics 
studies. 

In many aspects, the physics simulation included in Herwig++ is already superior to that in 
the FORTRAN HERWIG and our intention is that with the features just listed, the next major 
version release of Herwig++ will replace HERWIG as the recommended product for simulating 
hadron emission reactions with interfering gluons. 
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A Repository Commands 



The composition of the Repository is controlled through a simple command language, which can 
be used either interactively after calling Herwig++ read without any arguments, or through 
input files, which can be provided as arguments to the Herwig++ read command. The following 
overview only describes the most important repository commands. Examples of input files using 
this command language can be found in the 

HERWIGPATH/share/Herwig++ 
HERWIGPATH/share/Herwig++/defaults 

directories. Please note that the repository allows for an internal filesystem-like structure of 
directories and entries. This does not, however, correspond to any physical files on the operating 
system. 

We first give the commands that affect the overall state of the Repository, followed by com- 
mands for navigating the filesystem-like structure, event generation, creating and modifying 
objects in the Repository, and finally some miscellaneous commands. We conclude with a brief 
example of using the filesystem-like structure of the Repository to obtain the parameter values 
used in a run. 

Repository state 

save file 

Save the current repository state. 

load file 

Load a repository. Replaces the current state, 
read file 

Read in additional commands from file. 
library lib 

Load the dynamic shared library lib immediately, making all classes in the library available. 
Repository tree 

All operations in this section affect the repository tree only, not the file system, 
pwd 

Print the current directory path, 
cd dir 

Change the current directory to dir. 

mkdir dir 

Make a directory called dir. 

Is [dir] 

List the entries in the current directory or in dir. 
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rmdir dir 

Remove an empty directory. 



rrmdir dir 

Remove a directory and all its contents recursively. 

Event generation 

run run-name generator 

Run the generator object for the pre-set number of events. Files are saved under the label 
run-name. 

saverun run-name generator 

Save a generator as a file run-name. run, ready to use with Herwig++ run. 

Classes, objects, interfaces 

create classname name [library] 

Create a new object of C++ class classname and store it under name. Optionally, specify 
the name of the library file containing the class. 

mv old-name name 

Rename a repository object. 

cp old-name name 

Copy a repository object. The copy's interfaces will be identical to the original's at the 
time of copying, but can then be set independently. 

rm name 

Remove name from the repository. 

get interface 

Get the current value of an interface. 

set interface value 

Set the value of an interface. This can be either a numerical value, the name of an object 
in the Repository, or a defined key word for a Switch, set can also be used to set the value 
of a member of an interface vector. 

insert vector-interface [index] value 

Insert a value into a vector of interface parameters. 

erase vector-interface[index] 

Remove a value from a vector of interface parameters. 

describe object [interface] 

Describes object and lists its interfaces, or describes an interface. 
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Miscellaneous commands 

setup object args. . . 

Passes args to object's own setup function 34 . 

decaymode tag BR active? decay er 

Register a decay mode where tag is a semicolon-delimited description of a decay, using the 
repository particle names, such as piO->gamma,e-,e+;, BR is the mode's branching ratio, 
active? is either 1 or 0, indicating whether this decay mode is active or not, and decay er 
is the object that handles the generation of the kinematics for this decay mode. 

makeanti particlel particle2 

Register particlel and particle2 to be a particle-antiparticle pair. 

def aultparticle particle [particle . . .] 

Register particles as default particles, only these particles are used with every event gener- 
ator. 

A.l Example 

This is a brief example of using the Repository to extract the values of the default kinematic cuts 
on particles produced in the hard scattering process. Many more complicated tasks can also be 
performed. 

While we expect that the most common way of using the Repository will be changing the . in 
file for the relevant collider it is sometimes useful to browse the directory-like structure to check 
the parameters being used. 

The filesystem-like structure of the Repository can be explored using 

Herwig++ read 

which gives access to a command-line prompt. The current directory will be the last one used in 
the default Herwig++ Repository, currently /Herwig/ Analysis. Typing Is will give a list of the 
AnalysisHandler objects that have been created to analyse events generated by Herwig-I— h 

The objects that supply the kinematic cuts are in the directory /Herwig/Cuts and can be 
listed using 

cd /Herwig/Cuts 
Is 

which will list the following objects 

EECuts 

JetKtCut 

LeptonKtCut 

Mass Cut 

PhotonKtCut 

QCDCuts 

TopKtCut 

34 Used e.g. for particle data as 
setup particle ID PDGname mass width cut ctau charge colour spin stable 
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The QCDCuts and EECuts objects are the main objects that impose the cuts for hadron-hadron 
and lepton-lepton events respectively. Repository commands can now be used to get information 
about the objects and their parameters, for example 



describe QCDCuts 
describe QCDCuts : OneCuts 
get QCDCuts: OneCuts 

will give a brief description of the QCDCuts object and its interfaces, followed by the description 
of the OneCuts interface and the list of objects used to give the cuts on individual particles, or 
groups of particles. 

The JetKtCut object is used to impose cuts on partons (the quarks other than the top quark, 
and the gluon). The value of the cut on the the transverse momentum of the partons can be 
accessed and increased from the default value of 20 GeV to 30 GeV using 

get JetKtCut :MinKT 

set JetKtCut :MinKT 30.*GeV 

A new event generator file with this changed cut could now be written to file using 

saverun LHCnew LHCGenerator 

for the LHC. 
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B Examples 



This appendix contains a number of examples of using Herwig++. Example input files for 
Herwig++ are also supplied in the directory 

HERWIGPATH/share/Herwig++/ 

where HERWIGPATH is the location of the Herwig++ installation. There are examples for e + e~ 
collisions at LEP and ILC energies and hadron-hadron collisions at the Tevatron and LHC, as 
well as examples of using the different BSM models included in Herwig++. 
These can all be run with 

Herwig++ read Collider. in 

Herwig++ run -N no_of .events Collider. run 

where Collider, in is one of the example input files. The first read stage reads the input file 
and persistently writes the EventGenerator object it creates into the Collider. run file for future 
use. The second run stage then uses this persistently stored generator to generate no_of .events 
events. 

The default parameters for the generator have already been pre-set using the files contained 
in the directory 

HERWIGPATH/share/Herwig++/defaults 

and used to build the HerwigDefaults.rpo Repository file distributed with the release. Most users 
will not need to rebuild this file, but may need to look at the default parameters contained in 
the files used to build it. 

More information on running Herwig++ can be found on the wiki and in Appendix A. 

The remainder of this appendix is designed to illustrate how these input files can be adapted 
to simulate the physics scenario of interest to the user by changing the hard processes, cuts, etc. 
All of the examples, together with the source code, can be obtained from our wiki, where new 
examples will also be added in the future. Several of the examples assume that hadron-hadron 
collisions are being generated. If you are simulating lepton-lepton collisions replace LHCGenerator 
with LEPGenerator. In addition a number of useful examples can now be found in the Contrib 
directory of the release. 

B.l Switching parts of the simulation off 

In some cases it may be useful to switch off certain stages of the simulation. The most simple 
way to do that is by assigning NULL pointers to the appropriate StepHandlers of the EventHandler. 
The following statements have to be added to the Generator . in file used. 

cd /Herwig/EventHandlers 

set LHCHandler : CascadeHandler NULL 

set LHCHandler :HadronizationHandler NULL 

set LHCHandler :DecayHandler NULL 
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to switch off the parton shower, hadronization and hadronic decays. For e + e~ collisions the 
corresponding EventHandler is called LEPHandler. In e + e~ collisions it is possible, although not 
recommended, to switch the shower off while still hadronizing the event. This is not possible in 
hadron collisions because the decay of the hadronic remnant, which must occur before the event 
can be hadronized, is currently handled by the shower module. 

The Shower step can be controlled in more detail: Initial-state radiation can be turned off 
using 

set /Herwig/Shower/SplittingGenerator : ISR No 
Final-state radiation can be turned off using 

set /Herwig/Shower/SplittingGenerator : FSR No 
Multiple interactions can be turned off using 

set /Herwig/Shower/ShowerHandler :MPIHandler NULL 

By default Herwig++ now uses a multiple scattering model of the underlying event. If you 
wish to use the old UA5 model, which we do not recommend for realistic physics studies, you 
should first turn off the multiple scattering model and then enable the UA5 model 35 : 

set /Herwig/Shower/ShowerHandler :MPI No 
cd /Herwig/Hadronization/ 

set ClusterHadHandler :UnderlyingEventHandler UA5 
B.2 Setup for minimum bias runs 

With the introduction of the new underlying event model from Herwig++ 2.3 onwards, we are 
able to simulate minimum bias events. However, the default setup is tailored for simulating 
underlying event activity in hard scattering events. We therefore have to change the settings in 
a few places. First, the specific minimum bias matrix element, MEMinBias, has to be selected 

cd /Herwig/MatrixElements 

insert SimpleQCD:MatrixElements [0] MEMinBias 

Then, the cuts have to be altered to reflect the fact that we are using virtually any trigger 
selection. 

################################################## 

# Cuts on the hard process. MUST be ZERO for MinBias 
################################################## 

cd /Herwig/Cuts 

set JetKtCut:MinKT 0.0*GeV 

set QCDCuts:MHatMin 0.0*GeV 

# Need this cut only for min bias 
set QCDCuts:XlMin 0.01 

set QCDCuts:X2Min 0.01 

35 It should be remembered that there is a difference between the name of the class used to create objects in the 
Repository and the names of the objects, here ClusterHadHandler is the name of the ClusterHadronizationHandler 
object used by default in Herwig++ to perform the hadronization. 
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The remaining settings that need to be changed are the ones in the underlying event model, i.e. 
in the MPIHandler. The parameter that tells the model which hard process is identical to QCD 
jet production has to be set to zero as our primary hard process is QCD-like, 

set /Herwig/UnderlyingEvent/MPIHandler : IdenticalToUE 

The settings so far enable the new underlying event model with both perturbative and non- 
perturbative contributions. For completeness we report the switches that are available to turn 
the soft model off. In this case the simulation is identical to any version 2.x before 2.3. 

set /Herwig/UnderlyingEvent/MPIHandler : sof tint Yes 

The modification that calculates the overlap function of soft scatters from the elastic t-slope data 
can be controlled by 

set /Herwig/UnderlyingEvent/MPIHandler : twoComp Yes 

B.3 Simulation of several hard processes in one event 

In this section we show how several hard processes can be simulated in Herwig++. To achieve 
that, the hard processes have to be specified along with the cuts that should be used for them. 
We will choose the example of like-sign W production as illustrative example. 
We start with the conventional hard process and its cuts, where we select only W + production 
and decay to e + z/ e 

cd /Herwig/MatrixElements 

insert SimpleQCD:MatrixElements [0] MEqq2W2ff 

# only W+ 

set MEqq2W2ff :Wcharge 1 

# only positrons 

set MEqq2W2ff : Process 3 

To modify the cuts on that process we have to change the following 

cd /Herwig/Cuts 

set LeptonKtCut:MinKT 5.0*GeV 

set LeptonKtCut :MaxEta 2.5 

# inv mass cut on lepton pairs 
set MassCut:MinM 0.*GeV 

Now we can start adding additional hard processes. This is done in the MPIHandler class and 
always needs a SubProcessHandler that has a MatrixElement assigned and a compulsory reference 
to a Cuts object. This reference can be an existing one, in the case where we want to use identical 
cuts for the processes, but can also be an independent instance. To create such an independent 
instance we do 
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cd /Herwig/UnderlyingEvent/ 

# cut on pt . Without a specific matcher object, it works on all particles 
create ThePEG : : SimpleKTCut DPKtCut SimpleKTCut . so 

set DPKtCut :MinKT 10 
set DPKtCut : MaxEta 2.5 

# create the cuts object for DPI 
create ThePEG: : Cuts DPICuts 

# This should always be 2*MinKT! ! 
set DPICuts :MHatMin 20 

insert DPICuts : OneCuts DPKtCut 

We first created an instance of the class SimpleKTCut to implement a cut on the transverse 
momentum and rapidity of the outgoing particles. This instance then has to be assigned to the 
instance of the Cuts object, DPICuts. To create a valid SubProcessHandler we have to specify 

cd /Herwig/UnderlyingEvent/ 

create ThePEG :: SubProcessHandler DPI 

insert DPI :MatrixElements /Herwig/MatrixElements/MEqq2W2f f 
set DPI : PartonExtractor /Herwig/Partons/QCDExtractor 

We have assigned the reference to the same MatrixElement instance and therefore also have W + 
production and decay to positrons. The PartonExtractor is needed to extract the partons from 
the beam particles but is always the reference to the QCDExtractor. 

The last step now is to assign the subprocess and cuts instance to the MPIHandler and select the 
multiplicity at which they should appear. In our case this is of course simply one, but in the case 
of e.g. 6-quark pairs or QCD dijets this may be different from one. 

cd /Herwig/UnderlyingEvent/ 

# set the subprocesses and corresponding cuts 

# is reserved for the "usual" underlying events 

# Each SubProcessHandler must be accompanied by a Cuts object! 
insert MPIHandler : SubProcessHandlers 1 DPI 

insert MPIHandler : Cuts 1 DPICuts 

# now set what multiplicities you want, index = means the first 

# ADDITIONAL hard process 

# this is SubProcessHandler 1 with multiplicity 1 
insert MPIHandler : additionalMultiplicities 1 

B.4 Changing particle properties 

In Herwig++ each particle's properties are contained in a ParticleData object. This has a num- 
ber of interfaces that can be used to change the properties. The files leptons . in, quarks . in, 
bosons . in, mesons . in, baryons . in and diquarks . in, which can be found in the 
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HERWIGPATH/share/Herwig++/def aults directory, set up the default properties of each par- 
ticle type. The names of the ParticleData objects in the Repository can be found in these input 
files or by browsing the /Herwig/Particles directory in the Repository using Herwig++ read. 

All properties can be changed in the input file for an event generator. For example to change 
the mass of the top quark to 170 GeV the following lines should be added 

set /Herwig/Particles/t : NominalMass 170.*GeV 

By default, the properties of particles and their antiparticles are forced to be the same so this 
will change the mass of both the top and antitop. 
The neutral pion can be set stable using 

set /Herwig/Particles/piO : Stable Stable 
B.5 Changing some simple cuts 

In many cases it will be important to specify particular cuts on the hard process. The default 
values for all cuts in Herwig++ are given in the file 36 Cuts . in. Here we give a number of examples 
of changing the cuts. 

For example, in order to change the minimum k± for a parton produced in the hard process 
to 30 GeV one should add 

set /Herwig/Cuts/JetKtCut:MinKT 30.0*GeV 

The pseudorapidity cut on hard photons can be changed to \rj\ < 4 with 

set /Herwig/Cuts/PhotonKtCut:MinEta -4. 
set /Herwig/Cuts/PhotonKtCut:MaxEta 4. 

and the cut on the minimum invariant mass of the hard process can be increased to 50 GeV with 

set /Herwig/Cuts/QCDCuts:MHatMin 50.*GeV 

If one wants to restrict the invariant mass of the final state in lepton pair production, however, 
one should use the class V2LeptonsCut, our default instance of this is called MassCut. In this case 
one has to specify 

set /Herwig/Cuts/MassCut:MinM 20.*GeV 
B.6 Setting up an AnalysisHandler 

Creating a new AnalysisHandler requires the following steps, which should be done in a new 
directory outside the Herwig source tree: 

1. Create skeleton class files. This can be done in emacs by loading a Lisp script that can be 
found at THEPEG_PREFIX_PATH/share/ThePEG/ThePEG . el. 

36 This can be found in the directory HERWIGPATH/share/Herwig++/def aults 
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2. Now invoking M-x ThePEG-AnalysisHandler-class-f iles queries the user for some input 
and interactively creates the necessary files for an AnalysisHandler. These are the questions 
asked: 

(a) Class name: 

Use for example MyName : :Foo. It is useful to use a namespace (replacing MyName with 
your name, of course) 

(b) Base class name: 

The right answer is already suggested: AnalysisHandler 

(c) include file for the base class: 

Also filled out already 

(d) Will this class be persistent (y or n) 

If persistent members are needed: y otherwise n. n is appropriate here. 

(e) Will this class be concrete (y or n) 

The answer y is appropriate unless you're writing an abstract base class. 

This will create the following files: 
Foo.h, Foo.cc. Save all buffers now. 

3. If actions need to be performed as part of the initialization (e.g. booking histograms) or 
termination (e.g. writing results to disk), the required class methods can be automatically 
created by the same Lisp script: 

(a) First the declaration of the methods. Go to Foo.h where it says 

// If needed, insert declarations of virtual function defined in the 
// Interf acedBase class here (using ThePEG-interf aced-decl in Emacs) . 

and in emacs use M-x ThePEG-interf aced-decl. This will insert the declaration of 
the methods needed. 

(b) To insert the implementation of these methods, go to Foo.cc where it says 

// If needed, insert default implementations of virtual function defined 

// in the Interf acedBase class here (using ThePEG-interf aced-impl in Emacs) . 

and start M-x ThePEG-interf aced-impl. 

4. There is one important check left. Every class that can be administered by ThePEG has to 
specify a static function returning the name of the library that the class is stored in. This 
has to agree with the library name in the Makefile. In our case it is: 

static string library () { return "Foo.so"; } 

By default it is set to the name of the class, i.e. Foo . so in our case, but may need changing 
if you are linking several classes into one library. 

5. A fully working AnalysisHandler, which currently has no functionality, is now implemented. 
A Makefile to compile it is supplied with the release. Copy it to your working directory 
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cp HERWIGPATH/share/Herwig++/Makef ile-UserModules Makefile 

It will create a shared library object named after the . cc filename, e.g. Foo . so. 

6. The class can now be compiled by invoking make. This command should terminate suc- 
cessfully. 

7. Calling the newly created class requires copying an appropriate Generator, in file into 
your directory from HERWIGPATH/share/Herwig++ and modifying it with the following 

statements 

cd /Herwig/Analysis 

create MyName::Foo foo Foo. so 

insert /Herwig/Generators/LHCGenerator : AnalysisHandlers foo 

which will create an instance of the new class Foo and then insert it at position in 
the vector of references to AnalysisHandlers. It is always safest to insert the newly created 
AnalysisHandler as the first entry in the list unless you are sure of how many AnalysisHandlers 
have already been inserted. 

B.7 Usage of ROOT 

In the following we will show two examples of an AnalysisHandler that will use ROOT [167] output. 
Please refer to Appendix B.6 for the generic instructions on setting up an analysis. Here, we will 
only mention specific code snippets, which should be inserted in the appropriate locations. 
The short description of what has to be done is: 

1. create a new class derived from AnalysisHandler; 

2. implement the functionality required; 

3. compile a library from it; 

4. create a Generator, in file where this AnalysisHandler is called and run it. 

Points 1 and 3 are universal for every AnalysisHandler and are described in Appendix B.6. How- 
ever, the corresponding library and include statements for ROOT have to be added: First copy 
the Makefile 

cp HERWIGPATH/share/Herwig++/Makef ile-UserModules Makefile 

and then add the following lines 

ROOTCFLAGS := $ (shell root-config — cflags) 

ROOTGLIBS := $ (shell root-config — glibs) 

ROOT = $ (ROOTCFLAGS) $ (ROOTGLIBS) 

Finally the line containing the compilation command has to be changed to include the content 
of the ROOT variable: 
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°/ .so : %.cc °/„.h 

$(CXX) -fPIC $(CPPFLAGS) $( INCLUDE) $(ROOT) \ 
$(CXXFLAGS) -shared $< -o $@ 

A shared library with your code will be created in the directory where you execute make. You 
need to make sure that the ROOT libraries can be found at run-time. On Linux systems you can 
add paths to the libraries to the environment variable $LD_LIBRARY_PATH. 

B.7.1 Root histograms 

The goal of this example is to write an AnalysisHandler that writes the charged particle multiplicity 
per event to a histogram and saves it as an encapsulated postscript (eps) file. This is only a 
minimal example of the use of ROOT in the analysis of Herwig++ events. It may for example 
be more useful to write the histogram to a file, but we leave this to the user as it is beyond the 
scope of this manual. 

First a new AnalysisHandler has to be created, as described in Appendix B.6. After setting 
up the necessary files, the new functionality can be implemented: 

• Foo.h 

In the header file, several additional include files have to be specified 

#include "ThePEG/EventRecord/Event .h" 
#include "ThePEG/EventRecord/Particle . h" 
#include "ThePEG/EventRecord/StandardSelectors .h" 

#include "THIF.h" 
#include "TCanvas.h" 

The first one is needed to access methods of the event class. The next two are needed 
because we will use the Particle class to check for particle properties. The last two ones 
are ROOT headers of histograms and a canvas to draw the histogram on. 

The histogram should be available as a member of this new class, because information on 
every event has to be stored in it. A pointer to the histogram as private member variable 
of the class can be used for that purpose: 

private : 
/** 

* A pointer to a Root histogram 
*/ 

TH1F* histo; 

• Foo.cc 

The histogram should be booked in 

void Foo : : doinitrunO 
with the following commands: 
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histo = new THlF("test", "charged multiplicity", 150, 0, 600); 
histo->SetXTitle("N_{ch}") ; 
histo->SetYTitle(" events") ; 

In 

void Foo : : dof inish() 

the histogram is drawn on a canvas and saved to disk. Finally the pointers are freed: 

TCanvas *can = new TCanvas("plot" , ""); 
histo->Draw() ; 
can->SaveAs ( "plot . eps " ) ; 
delete can; 
delete histo; 

All that remains is the actual filling of the histogram. This functionality will be added to 
the method 

void Foo :: analyze (tEventPtr event, long, int loop, int state) { 
if ( loop > M state != II ! event ) return; 
/** create local variable to store the multiplicity */ 
int mult (0) ; 

/** get the final-state particles */ 
tPVector particles=event->getFinalState() ; 
/** loop over all particles */ 

for (tPVector :: const_iterator pit = particles .begin() ; 
pit != particles . end() ; ++pit){ 
/** Select only the charged particles */ 
if( ChargedSelector :: Check (**pit) ) 
++mult ; 

} 

histo->Fill(mult) ; 

} 

The test in the first line is recommended for all simple AnalysisHandlers. The meaning of 
loop and state can be obtained from the Doxygen documentation of the AnalysisHandler 
class. 



B.7.2 rtuple with TTree 

If you are working with ROOT already, you can store events in an rtuple directly. This example 
shows how to define an AnalysisHandler that prepares an rtuple with ROOT TTree. It is extracted 
from a more detailed example, available from the wiki, for analysing four-6 events at LEP. 

• Foo.h 

First, add the needed ROOT header files to your header file for declaration of all ROOT 
classes you are going to use. In this case: 
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#include "TTree.h" 
#include "TFile.h" 

Add TTree and TFile objects to the private part of the class: 

private : 
// ROOT Tree 
TTree * theTree; 
// ROOT File 
TFile * theFile; 

Define all the variables and arrays that will be kept in the ROOT tree: 
private : 

// ROOT tree internal arrays and variables 
int Nentry, Nqurk, Nhdrn; 
int Kf [16] , Kp[16] ; 
double Wgt , Alphas ; 
double Qscl [4] ; 

double Px[16], Py[16] , Pz[16], P0[16]; 
• Foo.icc 

Methods for TTree booking and the writing of the TFile to disk should be called in 
doinitrunO and dofinishO respectively. Add the following lines to doinitrun(): 

LEPbbbbComparison: :doinitrun () { 

// create ROOT Tree 

theTree = new TTree ("bbbb" , "my Analysis root tree", 1); 
if (! theTree) { 

cerr << "ROOT tree has not been created. . .\n" ; 

return; 

} 

// create ROOT File 

theFile = new TFile (outname, "RECREATE") ; 
if (! theFile) { 

cerr << "ROOT file has not been created. . .\n" ; 

return; 

} 

theTree->SetDirectory (theFile) ; 
// define ROOT Tree branches/leaves 
theTree->Branch ("Nentry", feNentry, "Nentry/I"); 
theTree->Branch ("Nqurk", feNqurk, "Nqurk/I"); 

theTree->Branch ("Pz", Pz, "Pz [Nentry] /D") ; 
theTree->Branch ("PO", PO, "PO [Nentry] /D") ; 
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} 



The last parameter in each command theTree->Branch() should be equal to "Name/Type" 
of each variable, e.g. I — > int, D — > double, etc. (Information on other types can be found in 
the ROOT manual). Final commands should be placed in LEPbbbbComparison: : dof inishO 
So, add the following lines to dofinishO: 

LEPbbbbComparison: : dofinishO { 

theTree->GetCurrentFile () ; 
theTree->Write() ; 
theFile->Close() ; 

cout << "ROOT file has been written on disk" << endl; 
} 

After that, the class will keep theTree in theFile and write theFile to disk. 
• Foo.cc 

All the TTree variables should be set in analyze (...). As soon as all the variables have 
the right values for analysing an event, execute the FillO method for theTree. 

void Foo :: analyze (tEventPtr event, long, int loop, int state) { 

// Fill TTree record 

if (2 < bquark.size ()) { 
theTree->Fill() ; 

} 

} 

B.8 Using BSM models 

There are example files installed in HERWIGPATH/share/Herwig++ that show how to use the 
implemented BSM physics modules. Each one is labelled Generator-Model . in. Also associated 
with each BSM physics module is a .model file that is required to run with a specific module but 
otherwise does not need to be touched by the user. The easiest method to run a BSM physics 
module is to copy the Generator-Model . in file that is appropriate to the collider and model 
under study and make the necessary changes there. 

B.8.1 MSSM 

To generate a process in the MSSM, first decide on the accelerator to use, the LHC for example, 
and then copy MSSM. model and LHC-MSSM. in files to the location where Herwig++ will be used. 
LHC-MSSM . in contains the settings that a user can manipulate, the default settings are for squark 
production at the LHC. To change this to gluino production one should delete the lines 
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insert HPConstructor : Outgoing /Herwig/Particles/~u_L 
insert HPConstructor : Outgoing 1 /Herwig/Particles/~u_Lbar 
insert HPConstructor : Outgoing 2 /Herwig/Particles/~d_L 
insert HPConstructor : Outgoing 3 /Herwig/Particles/~d_Lbar 

and insert the line 

insert HPConstructor : Outgoing /Herwig/Particles/~g 

A SUSY model requires a spectrum file to set the masses and couplings. This file is produced 
using a spectrum generator 37 . The name of the file, e.g. spectrum, spc, is set via the command 

setup MSSM/Model spectrum. spc 

If the decay table is in a separate file to the spectrum then a second setup line should be used 
to supply this file name. 

The next step is to set up the particles that will require spin correlations included in their 
decays. This is achieved through the DecayParticles interface. In the example of gluino 
production firstly one should remove the lines 

insert NewModel : DecayParticles /Herwig/Particles/~d_L 
insert NewModel : DecayParticles 1 /Herwig/Particles/~u_L 
insert NewModel : DecayParticles 2 /Herwig/Particles/~e_R _ 
insert NewModel : DecayParticles 3 /Herwig/Particles/~mu_R- 
insert NewModel : DecayParticles 4 /Herwig/Particles/~chi_10 
insert NewModel : DecayParticles 5 /Herwig/Particles/~chi_20 
insert NewModel : DecayParticles 6 /Herwig/Particles/~chi_2+ 

and then insert the line 

insert NewModel : DecayParticles /Herwig/Particles/~g 

This will generate spin correlations in the decay of the gluino but not in the subsequent decays 
of its children. Assuming these too are required then additional lines containing all of unstable 
products in the cascade decays are needed. 

insert NewModel : DecayParticles /Herwig/Particles/~g 
insert NewModel : DecayParticles 1 /Herwig/Particles/~d_L 
insert NewModel : DecayParticles 2 /Herwig/Particles/~u_L 

The rest of the settings in the file deal with general parameters for the run. Herwig++ can 
then be run as described at the beginning of this appendix. 



37 Some of these are listed at http : / /home . f nal . gov/ ~skands/slha/ 
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B.8.2 MUED 



The MUED model works in a similar fashion to the MSSM but with some important differences 
due to the unavailability of spectrum and decay generators. The mass spectrum is calculated by 
Herwig++ once the main parameters have been set via the interfaces 

set MUED/Model : InverseRadius 500.*GeV 
set MUED/Model : LambdaR 20 

and optionally 

set MUED/Model :HiggsBoundaryMass 0.*GeV 

Similarly to above the file LHC-MUED . in should be copied to a new file i. e. mymued . in and the 
relevant parameters changed there. 

The specification of the hard process is done in the same manner as above using the particle 
content of the MUED model. As there are no decay table generators for UED the possible pertur- 
bative decays are calculated automatically for the particles specified through the DecayParticles 
interface. It is advisable to leave the list as it stands in the file as then all of the necessary decays 
modes for the parents that are children in cascade decays will be created properly. 

Finally, the methods for running the generator are the same as above. 



B.8.3 RS Model 

Currently there are no matrix elements for the hard scattering that have tensor particles as 
external particles, they are only included as intermediates. The graviton can therefore only 
be included as a resonance. There is a special class designed to handle this as described in 
Appendix 5. 

The set up differs only slightly from the MSSM and MUED models. Using the example in 
LHC-RS . in, upon copying this to a new file, the lines 

insert ResConstructor : Incoming /Herwig/Particles/g 

insert ResConstructor : Incoming 1 /Herwig/Particles/u 

insert ResConstructor : Incoming 2 /Herwig/Particles/ubar 

insert ResConstructor : Incoming 3 /Herwig/Particles/d 

insert ResConstructor : Incoming 4 /Herwig/Particles/dbar 



insert ResConstructor : Intermediates /Herwig/Particles/Graviton 

insert ResConstructor : Outgoing /Herwig/Particles/e+ 
insert ResConstructor : Outgoing 1 /Herwig/Particles/W+ 
insert ResConstructor : Outgoing 2 /Herwig/Particles/ZO 
insert ResConstructor : Outgoing 3 /Herwig/Particles/gamma 

can be changed to suit the user's needs. The only parameter in this model is the cutoff scale and 
it is changed through the line 

set RS/Model:Lambda_pi 10000*GeV 

Again, running the generator follows the same steps as before. 
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B.8.4 Disabling Selected Decay Modes 

The decay modes for the new physics models do not exist prior to the Herwig++ read step so 
they cannot be simply commented out of a file. We therefore provide a separate mechanism to 
disable decay modes selected by a user. It is universal in that if a decay table is supplied with a 
SUSY model and a decay mode from this is specified, it will also be disabled. 

Using this mechanism simply requires adding information to the . in file. As an example we 
take the MUED model where we wish to disable the decay modes u\ — > 71, u and Z\ — > e°~, e + . 
In the relevant . in file, the lines 

cd /Herwig/NewPhysics/ 

insert DecayConstructor :DisableModes KKl_u_L->KKl_gamma,u; 
insert DecayConstructor :DisableModes 1 KKl_ZO->KKl_e_R-,e+; 

should be added, where the ordering of the decay products does not matter. The two characters 
following the parent particle are a dash (-) and a greater-than symbol (>). It should be noted 
that only the exact mode specified is disabled, i.e. in the above example Z\ — > ej + ,e~ would 
still be active. 



B.9 Intrinsic pt 

An example of a particular choice for the implementation of the intrinsic p? can be found in the 
default file Shower . in. 

set Evolver : IntrinsicPtGaussian 2.2*GeV 

As discussed in Appendix D, a Gaussian distribution for intrinsic p? has been implemented. 
The root mean square intrinsic pt of the Gaussian distribution required, a, is set using the 
IntrinsicPtGaussian parameter. The values for the intrinsic Pt are generated according to: 



d Pt — n exp 
ira z 



pry' 



a 



fB.r 



The default example above is for a Gaussian distribution with root mean square px of 2.2 GeV. In 
addition to this, there is the option of selecting an inverse quadratic distribution for the intrinsic 
Pt according to: 

vrln (1 + ^*P-\ 7 +p T 

where 7 is a constant and pr max is an upper-bound on the modulus of pt and makes the distri- 
bution normalizable. These parameters can be changed from their default values in Shower . in. 

set Evolver : IntrinsicPtGamma 0*GeV 
set Evolver : IntrinsicPtlptmax 0*GeV 

A mixture of both distributions can also be selected by setting a parameter (3 in Shower . in and 
is the proportion of the inverse quadratic distribution required and ranges between and 1. 

set Evolver : IntrinsicPtBeta 

Here the default setting is to generate the intrinsic px according to the Gaussian distribution 
only. 
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B.10 LesHouchesEventHandler 



In order to use partonic events generated by an external matrix element generator, a 
LesHouchesEventHandler object has to be created in the Repository. This object is supplied 
with at least one LesHouchesReader object. LesHouchesReader objects supply events in the Les 
Houches Accord (LHA) format [43] reading a file of events. 

Here we give an example of how to use LHA event files. The reading of the events is performed 
by the MadGraphReader class. This is not, however, limited to reading events generated by 
MadEvent [168] but can handle arbitrary event files in the Les Houches format. 

First, the libraries required must be loaded, 

library LesHouches . so 
library MadGraphReader . so 

Suppose the event file is called myEvents . lhe 38 . We will assume it contains some process of 
interest at the LHC First, a LesHouchesReader object needs to be created and given the name 
of the file: 

cd /Herwig/EventHandlers 

create ThePEG: : MadGraphReader myReader 

set myReader :FileName myEvents. lhe 

In principle, the information needed to generate full events, i.e. beam energies, incoming 
particles and parton distributions, is extracted from the event file, but may also be set explicitly. 
For these switches, see the interface documentation of the LesHouchesReader and MadGraphReader 
classes, respectively. 

In case files with unweighted events not generated by MadEvent are used, the LesHouchesReader 
needs to be assigned an event cache to gain information on the event sample. If, for example, 
events should be cached in a file named cacheevents . tmp the following setting should be used: 

set myReader : CacheFileName cacheevents . tmp 

The cuts on the hard process cannot, in general, be extracted from event files. If the interface 
value 

set myReader : InitCuts 

is assigned, the LesHouchesReader object expects to be given a Cuts object. For example, typical 
cuts for hadron collisions may be chosen: 

set myReader : Cuts /Herwig/Cuts/QCDCuts 

The use of cuts in Herwig++ is described in Appendix 3.5.3 and examples of changing them are 
given in Appendix B.5. If no Cuts object is assigned, the Cuts object assigned to the LesHouche- 
sEventHandler is used. 

Similar remarks apply to the use of parton distribution functions, which can be set explicitly 
using 

38 Thc LesHouchesReader class is also able to read in compressed event files, .lhe.gz. 
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set myReader : InitPDFs 

set myReader :PDFA firstBeamPDF 

set myReader :PDFB secondBeamPDF 

where firstBeamPDF and secondBeamPDF are PDFBase objects. Here, either the built-in PDF 
set or LHAPDF may be used, see Appendix B.ll. 

Next a LesHouchesEventHandler object has to be created. Objects of this class inherit from 
EventHandler and provide the same interfaces. The setup is therefore similar to the setup of 
a Standard EventHandler object, which needs to be equipped with showering, hadronization and 
decay handlers: 

create ThePEG: : LesHouchesEventHandler myLesHouchesHandler 

set myLesHouchesHandler : CascadeHandler /Herwig/Shower/ShowerHandler 

set myLesHouchesHandler :HadronizationHandler \ 

/Herwig/Hadronization/ClusterHadHandler 
set myLesHouchesHandler :DecayHandler /Herwig/Decays/DecayHandler 
set myLesHouchesHandler : PartonExtractor /Herwig/Partons/QCDExtractor 

A Cuts object that is applied to all processes may be set as for every EventHandler. Finally, 
the LesHouchesReaders from which the event handler should draw events have to be specified: 

insert myLesHouchesHandler : LesHouchesReaders myReader 
insert myLesHouchesHandler : LesHouchesReaders 1 myOtherReader 

An arbitrary number of readers may be used. 

A default or custom EventGenerator object can use the LesHouchesEventHandler object 
myLesHouchesHandler and a run file can be created from this event generator: 

cd /Herwig/Generators 

cp LHCGenerator myLesHouchesGenerator 

set myLesHouchesGenerator : EventHandler \ 

/Herwig/EventHandlers/myLesHouchesHandler 
saverun myLesHouches myLesHouchesGenerator 

The event generator can then be used as described at the beginning of Appendix B. 
B.ll Use of LHAPDF 

Herwig++ provides a built-in PDF set 39 . Other PDF sets may be used through the LHAPDF [170] 
interface of ThePEG. This section contains an outline of the use of LHAPDF. 
ThePEG has to be configured to use LHAPDF by adding the option 

— with-LHAPDF=/ path/to/LHAPDF/1 ib 

to the call of the configure script. Note that the full path to the LHAPDF libraries needs to 
be given. Once Herwig++ is built using ThePEG configured to use LHAPDF, PDF sets can be 
created easily in the Repository, for example the CTEQ6L set: 

39 The default PDF set in Herwig++ is the leading-order fit from the MRST'02 family [169]. 
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create ThePEG: :LHAPDF myPDFset 
set myPDFset :PDFName cteq61 . LHpdf 

set myPDFset :RemnantHandler /Herwig/Partons/HadronRemnants 
set /Herwig/Particles/p+:PDF myPDFset 
set /Herwig/Particles/pbar- :PDF myPDFset 

The custom PDF set myPDFset may also be used in a LesHouchesReader object, see Ap- 
pendix B.10. 

B.12 Use of a simple saturation model for PDFs 

A very simple modification that mimics parton saturation effects can be applied for any PDF by 
using the SatPDF class. The modification replaces xf(x) below xq by 



where XO and Exp are the changeable parameters. After copying an appropriate Collider, in 
to your local directory, adding the following lines before any run or saverun statement will enable 
the PDF modifications. 

################################################## 
# saturation modifications 

################################################## 
cd /Herwig/Partons 

create Herwig: : SatPDF SaturationMod HwSatPDF.so 
set SaturationMod : RemnantHandler HadronRemnants 

## Assign the pdf that should be modified: 

## use internal pdf 

set SaturationMod: PDF MRST 

## use lhapdf . This depends on the name you have 
## chosen for the LHAPDF set 
#set SaturationMod: PDF foo 

## may change XO: default is 1E-4 
#set SaturationMod :X0 1E-3 

## may change Exp: default is 
#set SaturationMod: Exp 1 

## Assign the modified pdf to the beam particles, 
## without this step the original pdf ; s will be used 
set /Herwig/Particles/p+ : PDF SaturationMod 
set /Herwig/Particles/pbar- : PDF SaturationMod 
cd /Herwig/Generators 
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C Contrib 



Starting with Herwig++ 2.3 we include a number of modules in the Contrib directory supplied 
with the release. In general code in this directory falls into one of four categories: 

1. Code, generally AnalysisHandlers, that was written by Herwig++ authors to test the imple- 
mentation of some new feature in Herwig++, which we expect that the vast majority of 
users will not need, however we now distribute it in case it is useful to some users; 

2. Code that depends on external libraries, for example Rivet, which we do not wish all of 
Herwig++ to depend on but may provide a useful interface for some users; 

3. Implementations of new models that are not sufficiently important to be included in the 
core Herwig++ code but may be of interest to some users; 

4. Code supplied by people who are not authors of Herwig++ that may be of use to some 
users. 

Currently the following modules are included: 
Analysis2 An alternative analysis framework; 

LeptonME Simple matrix elements for testing at the T(4s) resonance and testing Higgs decays; 
TauAnalysis Analysis of tau decays, used in Ref. [25]; 

DecayAnalysis Analysis of hadron decays, used to test the Herwig++ hadron decays; 
Rootlnterface An interface to the ROOT analysis framework; 
RivetAnalysis An interface to the Rivet analysis framework. 
RadiativeZPrime Implementation of the model of [171]; 

AnomalousHVV A simple implementation of anomalous Higgs boson couplings to vector 
bosons based on [172]; 

These modules are supplied with a simple Makefile, which will correctly link them against 
Herwig++ although if they interface to external packages some environment variables may have 
to be set by hand. 

In general these modules come with a much lower level of support and testing than the code 
in the main Herwig++ release, but if you do have problems, the authors of the contributed code 
will try to help where possible. Given the lower level of testing we request that you contact a 
module's authors before using its results in any publications. 
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D Tuning 

The default hadronization and shower parameters in Herwig++ have been tuned to a wide range 
of experimental data, primarily from LEP and 5-factory experiments. 

The following experimental data were used, with the exception of charm hadron spectra from 
the 5-factory experiments, all are from e + e~ experiments operating at the Z° peak: 

• the momentum spectra of charm hadrons, i.e. D* ±,Q , D ±,Q , Df, and A+, measured by the 
Belle collaboration away from the T(4S) resonance, [173]; 

• the momentum spectra of charm hadrons, i.e. D* ±)0 and D ±0 , measured by the CLEO 
collaboration away from the T(4S) resonance, [174]; 

• the weakly decaying 5-hadron fragmentation functions measured by the ALEPH [175] and 
SLD [176] collaborations; 

• four-jet angles measured by the ALEPH collaboration [177]; 

• the momentum spectrum of charged particles, charged pions, charged kaons and protons 
for all, light, charm and bottom quark events measured by the SLD collaboration [178]; 

• the momentum spectra for the production of 7r ± [179], K ± [179], p [179], A ++ [180], 
~*° [181], f 2 [182], / (980) [182], <\> [182], K*° [183], K° [184], vr [185], r) [185], r} 1 [185], 
p ± [185], u [185], 4 [185], a - [181], £*± [181], measured by the OPAL collaboration; 

• the multiplicity of charged particles measured by the OPAL collaboration [186]; 

• the momentum spectra for the production of p° [187] and D° [188] measured by the DELPHI 
collaboration; 

• the momentum spectrum of D*^ mesons measured by the ALEPH collaboration [189]; 

• the momentum spectrum of A baryons [190] and K*^ mesons [190] measured by the 
ALEPH collaboration; 

• the differential distributions y nm where an event changes from being an n to an m jet event 
according to the Durham jet algorithm, jet production rates and the average jet multiplicity 
as a function of the Durham jet measure measured by the OPAL collaboration [191]; 

• the differential jet rates with respect to the Durham jet measure measured by the DELPHI 
collaboration [192]; 

• the thrust, thrust major, thrust minor, sphericity, oblateness, planarity, aplanarity, C and D 
parameters, hemisphere masses, and jet broadening event shapes measured by the DELPHI 
collaboration [192]; 

• the rapidity, and transverse in and out of the event plane with respect to the thrust and 
sphericity axes measured by the DELPHI collaboration [192]; 
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• the average multiplicities of charged particles, photons, 7r°, p°, 7r + , p + , rj, u>, fi, K°, K*°, 
K*°, K+, K*+, rjf, 0, & D+, D*+, D°, D+, J/^, n°, p+, A++, E*~, A , S°, £+, £*+, 
S- fi- A+ f^f u ifj(2S), a+ taken from the PDG [49]; 

• the fractions of B°, B ± and 6-baryons from the Heavy Flavour Averaging Group (HFAG) [193]. 
The following parameters were tuned: 

1. the value of as at the Z° mass, AlphaMZ; 

2. the cutoff scale in the parton shower cutoffKinScale; 

3. the ConstituentMass of the gluon used in the hadronization model; 

4. the scale Qmin below which a non-perturbative treatment of as is used, the default is to 
set as to a constant below this scale; 

5. the maximum mass CIMaxLight above which clusters containing light quarks undergo 
cluster fission, see Eq. (7.1); 

6. the exponent CIPowLight controlling whether clusters containing light quarks undergo 
fission, see Eq. (7.1); 

7. the exponent PSplitLight controlling the masses of the daughter clusters for light quark 
clusters that undergo fission, see Eq. (7.2); 

8. the CISmrLight parameter, which controls the smearing of the direction of hadrons con- 
taining perturbatively produced light quarks, see Eq. (7.14); 

9. the weight PwtSquark for producing a strange quark- ant iquark pair in the hadronization; 

10. the weight PwtDIquark for producing a diquark-antidiquark pair in the hadronization; 

11. the relative weight SngWt for the production of singlet baryons; 

12. the relative weight DecWt for the production of decuplet baryons; 

13. the maximum mass CIMaxCharm above which clusters containing charm quarks undergo 
cluster fission, see Eq. (7.1); 

14. the exponent CIPowCharm controlling whether clusters containing charm quarks undergo 
fission, see Eq. (7.1); 

15. the exponent PSplitCharm controlling the masses of the daughter clusters for charm 
quark clusters that undergo fission, see Eq. (7.2); 

16. the CISmrCharm parameter, which controls the smearing of the direction of hadrons 
containing perturbatively produced charm quarks, see Eq. (7.14); 

17. the SingleHadronLimitCharm parameter, which controls the splitting of charm clusters 
to a single hadron above the threshold for producing two hadrons, see Eq. (7.15); 
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18. the maximum mass CIMaxBottom above which clusters containing bottom quarks un- 
dergo cluster fission, see Eq. (7.1); 

19. the exponent CIPowBottom controlling whether clusters containing bottom quarks un- 
dergo fission, see Eq. (7.1); 

20. the exponent PSplitBottom controlling the masses of the daughter clusters for bottom 
quark clusters that undergo fission, see Eq. (7.2); 

21. the CISmrBottom parameter, which controls the smearing of the direction of hadrons 
containing perturbatively produced bottom quarks, see Eq. (7.14); 

22. the SingleHadronLimitBottom parameter, which controls the splitting of bottom clus- 
ters to a single hadron above the threshold for producing two hadrons, see Eq. (7.15); 

The tuning was performed in a number of stages: 

• 200,000 events were generated at each of 2000 randomly selected parameter points for the 
first 7 parameters, which are sensitive to general properties of the events; 

• for the values of the first 7 parameters that gave the lowest \ 2 from the first scan 200,000 
events were generated for randomly selected values of parameters 8-11, which mainly con- 
trol the multiplicities of different hadron species; 

• for the values of the first 11 parameters that gave the lowest \ 2 from the second scan 
200,000 events were generated for randomly selected values of parameters 12-21, which 
mainly control the production of bottom and charm hadrons; 

• the parameters were then scanned about the minimum \ 2 point and the parameter that 
gave the largest reduction in the x 2 was adjusted to the value that gave the minimum value; 

• the scanning of parameters about the minimum was repeated until no significant improve- 
ment was found; 

• finally some parameters, particularly in the charm and bottom sector, that are not partic- 
ularly sensitive to the global \ 2 were adjusted to reduce the x 2 f° r observables sensitive to 
them. In practice the parameters 13-16 were adjusted to improve the quality of the fit to 
charm hadron multiplicities and spectra, the parameters 17-21 were adjusted to improve 
the quality of the fit to bottom hadron multiplicities and spectra, and the parameters 10-11 
were adjusted to improve the quality of the fit to baryon multiplicities and spectra. 

In each case 200,000 events were generated at both the Z° pole for the LEP observables and 
below the T(AS) resonance for the 5-factory observables. The \ 2 value included all the observ- 
ables but in order to increase the sensitivity to the particle multiplicities the \ 2 for the total 
particle multiplicities were multiplied by 10 when computing the global x 2 , and the total charged 
multiplicities at LEP by 100. 

The variation of the x 2 is shown in Fig. 8 for some of the parameters that are sensitive to the 
event shapes and production of hadrons containing the light, i.e. down, up and strange, quarks. 
The best fit point has a \ 2 — 6-4, with the increased weights for the hadron multiplicities and, 
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X 2 = 5.4 if all observables have unit weight. While this may seem too high a value, given the 
limited nature of the tuning it is not out of line with previous event generator tunings and the 
X 2 is about 4 times lower than before the tuning. 

In addition to the above, the option of including an intrinsic transverse momentum for partons 
within a hadron in hadron-hadron collisions has been implemented. It is chosen from the Gaussian 
distribution shown in Appendix B.9. For Drell Yan Z/W boson production at the Tevatron 
(y^ = 1.96 TeV), the best fit tune has an rms transverse momentum of 2.2 GeV [16]. For the 
CERN ISR experiment (y/s = 62 GeV) likewise, a best fit rms value of 0.9 GeV was obtained. 
Assuming a linear dependence of the rms value on ln(M/ ^/s) where M is the invariant mass, 
the corresponding value estimated for Z/W boson production at the LHC is within the range 
3.7 — 7.7 GeV. It is worth noting that the lower value of 3.7 GeV gives the best agreement with 
an alternative model [26], which introduces non-perturbative smearing during the perturbative 
evolution by modifying the implementation of ag- 
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Figure 8: Variation of the x 2 about the minimum points for the a) AlphaMZ, b) cutoffKinScale, 
c) CIMaxLight, d) PSplitLight, e) PwtSquark, and f) PwtDIquark parameters. The solid line 
shows the total x 2 i the dot-dashed line shows the x 2 f° r the particle multiplicities and the dotted line 
shows the x 2 f° r the event shape observables. In e) the dashed lines show the x 2 f° r observables sensitive 
to strange hadron production and in f ) the dashed lines show the x 2 f° r observables sensitive to baryon 
production. The vertical dashed lines show the final values of the parameters, described as default 
throughout this manual. In each figure, all other parameters are kept at their default values. 
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